Time  Resolution  of  Collapse  Events 
During  the  Propagation  of 
Ultraviolet  Light  Filaments 

THESIS 

Teresa  J.  Fondren,  First  Lieutenant,  USAF 
AFIT/GAP /ENP /08-M03 

DEPARTMENT  OF  THE  AIR  FORCE 
AIR  UNIVERSITY 

AIB,  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Patterson  Air  Force  Base,  Ohio 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the 
official  policy  or  position  of  the  United  States  Air  Force,  Department  of  Defense, 
or  the  United  States  Government. 


AFIT /GAP /ENP/08-M03 


Time  Resolution  of  Collapse  Events  During  the 
Propagation  of  Ultraviolet  Light  Filaments 


thesis 


Presented  to  the  Faculty 
Department  of  Engineering  Physics 
Graduate  School  of  Engineering  and  Management 
Air  Force  Institute  of  Technology 
Air  University 

Air  Education  and  Training  Command 
In  Partial  Fulfillment  of  the  Requirements  for  the 
Degree  of  Master  of  Science  in  Applied  Physics 


Teresa  J.  Fondren,  B.S.  in  Physics  and  in  Mathematics 
First  Lieutenant,  USAF 


March  2008 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


AFIT/GAP /ENP /08-M03 


Time  Resolution  of  Collapse  Events  During  the 
Propagation  of  Ultraviolet  Light  Filaments 


Teresa  J.  Fondren,  Bachelor  of  Science  in  Physics  and  in  Mathematics 

First  Lieutenant,  USAF 


Approved: 


Thomas  Niday  (Chairman 


//  Oct  07 

Date 


David  Weeks  (Member) 


//  o  7 

Date 


Matthew  BohM(Member) 


/  f  c?cA~ 

Date 


AFIT /GAP /ENP/08-M03 


Abstract 

Long  distance  propagation,  or  filamentation,  of  short,  intense  laser  pulses  is 
suggested  to  be  possible  through  the  balance  of  two  effects:  self-focusing,  when  a 
nonlinear  index  of  refraction  of  air  is  induced  by  high  intensities,  and  de-focusing, 
due  to  the  plasma  created  by  the  pulse.  Applications  for  filamentation  include 
areas  such  as  remote  sensing  and  directed  energy.  A  split-step  spectral  propaga¬ 
tion  simulation  is  used  to  model  the  behavior  of  a  high  intensity  ultraviolet  laser 
pulse  propagating  through  air.  Convergence  of  femtosecond  duration  collapses 
that  form  on  the  leading  edge  of  the  pulse  in  the  time  domain  is  achieved  with 
an  increase  in  the  multi-photon  ionization  coefficient.  Through  an  analysis  of  the 
relative  sizes  of  each  term  in  the  propagation  equation,  a  lack  of  plasma  present 
at  the  leading  edge  of  the  pulse  is  found  to  cause  these  collapses.  Results  for  a 
more  recent  value  of  the  electron-positive  ion  recombination  rate  are  compared 
to  results  from  a  higher  value  used  in  previous  works.  A  linear  stability  analysis 
shows  inherent  instability  of  the  pulses  in  all  cases.  The  inclusion  of  group  velocity 
dispersion  is  shown  to  increase  stability  at  high  temporal  frequencies  except  at 
zero  spatial  frequencies.  A  run  similar  to  an  experiment  that  produced  UV  fila¬ 
mentation  is  shown  to  propagate  70  cm  without  exceeding  numerical  limitations. 
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Time  Resolution  of  Collapse  Events  During  the 
Propagation  of  Ultraviolet  Light  Filaments 

I.  Introduction 

Motivation 

Dedicated  research  is  currently  being  done  on  the  long  distance  propaga¬ 
tion  of  intense  ultrashort  laser  pulses  through  the  atmosphere.  Applications  for 
this  include  remote  sensing  [1-3],  directed  energy  [4,5],  and  laser  induced  light¬ 
ning  [6,7].  All  of  these  applications  are  especially  interesting  to  the  military  for 
protecting  forces  from  harmful  environmental  effects,  shielding  them  from  enemy 
sensing  capabilities,  allowing  them  access  to  information  that  has  been  previously 
unattainable,  and  inflicting  damage  to  enemy  forces.  In  order  for  these  applica¬ 
tions  to  be  optimized,  two  things  are  required.  First,  to  further  the  effective  reach, 
the  range  of  propagation  through  atmospheric  phenomena  such  as  rain,  clouds, 
and  dust  particles  needs  to  increase.  Second,  to  be  most  effective,  more  power 
needs  to  be  transferred  with  the  pulses  all  the  way  to  the  target  [8]. 

Long  propagation  distances  of  laser  pulses  have  been  demonstrated  in  both 
the  infrared  (IR)  and  ultraviolet  (UV)  wavelength  ranges  due  to  hlamentation  of 
the  pulses  [9-11].  This  thesis  will  use  a  definition  of  hlamentation  or  filaments 
from  Couairon: 
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Filamentation  [refers  to]  a  dynamic  structure  with  an  intense  core  that 
is  able  to  propagate  over  extended  distances  much  larger  than  the 
typical  diffraction  length  while  keeping  a  narrow  beam  size  without 
the  help  of  any  external  guiding  mechanism  [12]. 

Another,  stricter  definition  is  also  used  in  the  literature  that  refers  to  filaments 
as  only  the  part  of  the  beam  propagation  where  the  pulse  generates  a  column  of 
weakly  ionized  plasma  in  its  wake  due  to  a  sufficiently  high  intensity  in  its  core. 

Filamentation  is  induced  when  the  power  of  a  pulsed  beam  is  above  a  cer¬ 
tain  critical  power  that  is  proportional  to  the  square  of  the  wavelength.  The  high 
power  initiates  a  nonlinear  effect  known  as  the  optical  Kerr  effect  which  causes 
self-focusing  of  the  laser  pulse.  The  Kerr  effect  introduces  an  intensity  dependent 
refractive  index  when  the  pulse  power  is  equal  to  or  greater  than  the  critical  power, 
creating  an  effective  focusing  lens  in  the  refractive  profile  of  the  air  [8].  While 
focusing,  the  pulse  eventually  reaches  a  high  intensity  as  it  overcomes  diffrac¬ 
tion.  This  causes  ionization  of  the  air  around  it,  forming  a  surrounding  plasma. 
The  plasma,  in  turn,  absorbs  photons  from  the  pulse  and  changes  the  index  of 
refraction  of  air,  resulting  in  a  decrease  in  pulse  intensity  and  an  effective  nega¬ 
tive  lens  that  defocuses  the  pulse.  Studies  have  suggested  that  balancing  these 
focusing-defocusing  effects  will  allow  propagation  over  long  distances  on  the  order 
of  several  Rayleigh  ranges  or  further  (up  to  kilometers  [13]),  possibly  even  through 
atmospheric  obscurants  [2,3,9-12,14], 

Much  of  the  research  involving  long-distance  propagation  of  light  filaments 
in  air  due  to  self-focusing  deals  with  wavelengths  in  the  IR  range  (>750  nrn).  A 


2 


main  attraction  to  this  wavelength  is  due  to  an  extra  effect  called  ultra-broadband 
radiation  that  light  filaments  of  these  wavelengths  can  generate.  This  radia¬ 
tion  emits  a  white-light  super-continuum  range  of  wavelengths  from  500  nm  to 
1200  nm,  which  can  be  used  to  remotely  monitor  trace  gases  and  aerosols  in  the 
air  [1,  15].  By  directing  such  a  continuum  of  wavelengths  at  a  target,  a  wide 
range  of  aerosols  could  be  detected  through  observing  the  absorption  of  certain 
wavelengths  from  the  continuum.  Early  warning  systems  for  bioterrorist  threats 
could  be  achieved  based  on  selective  remote  sensing  of  aerosols  with  the  ultra- 
broadband  radiation  induced  from  IR  hlamentation  since  aerosols  can  be  used  as 
carriers  for  toxic  species  to  be  absorbed  into  the  lungs.  Extending  the  Lidar  (light 
detection  and  ranging)  technique  with  nonlinear  and  multispectral  measurements, 
the  white-light  super-continuum  effect  could  allow  the  retrieval  of  the  abundance, 
particle  size,  and  refractive  index  of  a  diverse  set  of  aerosol  particles.  Previously, 
remote  identification  of  these  aerosols  was  extremely  difficult  due  to  the  variety 
of  their  size,  shape,  and  composition  [16]. 

Although  this  white  light  phenomenon  is  not  reported  for  UV  filaments  [9], 
using  wavelengths  of  10-400  nm  has  other  wanted  characteristics  for  hlamentation. 
Since  the  critical  power  needed  for  hlamentation  is  proportional  to  the  square  of 
the  wavelength,  UV  filaments  require  a  much  lower  peak  power,  with  an  on-axis 
intensity  on  the  order  of  1  TW/cnr,  in  comparison  to  100  TW/cm2  for  IR  [5,17]. 
This  results  in  three  possible  advantages: 


3 


•  Longer  UV  pulses  can  be  used  while  still  being  able  to  have  enough  power  for 
blamentation.  These  longer  pulses  should  be  able  to  travel  longer  distances 
since  the  filaments  can  contain  more  energy  [18]. 

•  Many  of  the  higher  order  nonlinear  effects  important  for  IR  ultra-short 
pulses  can  be  neglected  for  longer  UV  pulses. 

•  The  focus-defocus  balance  for  long  distance  propagation  could  become  solely 
intensity  dependent  allowing  the  increase  of  pulsewidth  to  increase  filament 
energies  [5]. 

These  potential  advantages  of  UV  pulses  create  a  demand  for  more  research  in¬ 
volving  UV  filaments.  There  is  a  need  to  better  understand  the  physics  behind 
UV  filament  propagation  in  order  to  exploit  its  usefulness. 

One  of  the  most  interesting  applications  of  this  hlamentation  is  in  the  area 
of  self-healing.  Work  has  been  done  to  investigate  the  effects  of  atmospheric 
obscurants  on  the  long  distance  propagation  of  laser  pulses.  Skupin  et  al.  [14] 
were  able  to  theoretically  show  that  single  filaments  can  survive  interaction  with 
obscurants  sized  up  to  two-thirds  that  of  the  filament  diameter.  This  agrees  with 
the  experimental  data  of  Courvoisier  et  al  [2],  The  durability  of  the  filaments 
is  due  to  the  ability  of  some  beam  components  that  are  kept  untouched  after 
collision  with  the  obscurant  to  refocus  and  continue  with  a  power  above  critical. 
Through  self-focusing  again  onto  the  beam  axis,  these  components  are  able  to 
replenish  the  filament  within  a  few  centimeters.  In  other  words,  instabilities  in 
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the  filament,  caused  by  the  high  powers  needed  to  overcome  the  critical  power 
required  for  filamentation,  allow  it  to  split  into  several  smaller  filaments  when  it 
encounters  an  obscurant.  Afterwards,  the  filaments  recombine  or  self-heal  and 
continue  along  the  original  path  with  only  a  10-15%  loss  in  power  [3].  This  could 
allow  long  range  sensing  to  be  done  through  thin  clouds  and  other  obscurants. 
This  ability  to  propagate  through  adverse  weather  would  be  invaluable  for  military 
remote  sensing  applications. 

Laser  induced  lightning  is  another  application  of  filamentation.  It  is  based 
on  the  plasma  wake  left  by  the  laser  as  the  filaments  propagate.  UV  laser  pulses, 
with  their  smaller  wavelength,  are  better  suited  for  this  application  since  they  can 
produce  a  longer  conducting  plasma  column  than  IR  laser  pulses  [12].  The  idea  is 
to  use  this  plasma  channel  created  by  the  filamentation  of  a  high-powered  short- 
pulsed  laser  to  trigger  and  guide  a  high-voltage  discharge  such  as  lightning  [6, 
16,19].  By  controlling  where  and  how  lightning  strikes  in  a  storm,  sensitive  sites 
such  as  electrical  installations  or  airports  could  be  protected  from  direct  strikes 
and  electromagnetic  perturbations.  Rodriguez  et  al.  were  able  to  demonstrate 
the  complete  guiding  of  a  discharge  produced  at  the  high-voltage  facility  of  the 
Technical  University  of  Berlin  [20]. 

This  same  idea  could  be  used  for  application  in  directed  energy.  Instead 
of  lightning  being  induced  and  guided,  the  plasma  channel  created  by  the  laser 
would  guide  a  purposed  electrical  discharge  directly  to  a  target  to  either  destroy 
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it  or  limit  its  capabilities.  In  reference  [4],  Davis  et  al.  suggest  a  four-step  process. 
First,  a  laser  creates  a  plasma  channel  that  has  low  conductivity  all  the  way  to 
the  target.  Second,  the  conductivity  of  the  channel  increases  by  several  orders  of 
magnitude  as  a  leader  is  propagated  down  it.  Third,  a  return  stroke  increases  the 
channel  conductivity  even  further,  allowing  swift  and  efficient  delivery  of  energy. 
Last,  the  ohmic  dissipation  along  the  line  drives  the  conductivity  up  even  further 
with  pulse  delivery.  This  four-step  process  could  allow  military  forces  to  send  a 
low  impedance  lightning-like  arc  strike  to  destroy  a  target. 

There  are  many  other  potential  applications  for  light  hlamentation  such  as 
generation  of  single  cycle  pulses  and  high  order  harmonics,  laser  induced  break¬ 
down  spectroscopy  (LIBS),  power  supply  for  high  speed  electric  vehicles,  and  even 
advances  in  propulsion.  A  detailed  summary  of  these  applications  and  the  research 
done  up  to  this  point  in  the  area  of  hlamentation  is  presented  by  Couairon  and 
Mysyrowicz  in  reference  [12]. 

Problem  Statement 

Although  there  has  been  a  thorough  amount  of  research  concerning  IR  light 
filaments  and  much  of  the  science  is  understood  and  can  be  reproduced  numer¬ 
ically  and  experimentally,  the  details  of  UV  hlamentation  are  still  debated.  In 
order  to  tap  into  the  benefits  of  using  UV  pulses,  a  better  understanding  of  what 
is  physically  happening  in  the  long  distance  propagation  of  UV  hlaments  is  needed. 
Many  studies  [5,8,9,11,21-23]  have  investigated  the  stability  of  UV  pulses  with 
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respect  to  the  balancing  of  self-focusing  and  diffraction  terms  in  the  propagation 
equation.  There  have  been  experimental  results  showing  propagation  of  UV  pulses 
to  distances  up  to  twelve  meters  [5].  The  question  remaining  and  being  discussed 
in  papers  on  the  subject  of  UV  hlamentation  is  how  these  filaments  are  able  to 
travel  these  distances  with  stability. 

Currently,  long  pulse  models  of  UV  hlamentation  exhibit  spatial  and  tem¬ 
poral  instabilities  within  the  pulse,  causing  it  to  breakup,  possibly  limiting  its 
propagation  distance  [21,23].  In  Niday’s  long  pulse  propagation  results,  the  pulses 
form  collapse  events  on  the  leading  edge  of  the  pulse.  Figure  1.1  is  taken  from  Ni¬ 
day’s  dissertation  [21]  and  shows  the  development  of  collapse  events  on  the  leading 
edge  of  the  pulse  and  the  instabilities  exhibited  as  it  propagates.  These  peaks  are 
three  magnitudes  of  order  shorter  in  time  than  the  original  pulse,  thus  making 
it  difficult  to  numerically  model  the  whole  pulse  while  having  enough  resolution 
to  be  able  to  understand  the  details  of  the  collapses.  Propagation  to  distances  of 
about  one  meter  with  this  code,  in  some  cases,  creates  an  overflow  of  the  numerical 
parameters  where  the  computational  area  (spatial  grid)  is  not  fine  enough.  In  the 
spatial  frequency  representation,  non-zero  values  start  to  exist  at  the  boundaries 
that  are  then  wrapped  around  to  the  opposite  edge  causing  non-physical  results 
that  lead  to  the  disintegration  of  the  pulse  starting  from  the  trailing  edge.  Higher 
computing  power  that  can  keep  track  of  a  larger  grid  or  better  coding  that  can 
allow  for  a  larger  grid  with  the  same  amount  of  computing  power  is  needed  to  in- 
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vestigate  the  causes  of  these  leading  edge  collapse  events  and  be  sure  to  eliminate 
effects  caused  by  numerical  grid  overflow.  Understanding  these  details  will  lead  to 
a  better  understanding  of  the  interaction  between  focusing  and  de-focusing  terms 
in  the  propagation  equation,  opening  the  door  to  applications  of  long  distance 
propagation  of  UV  pulses. 


Figure  1.1:  Propagation  of  a  finite  duration  pulse  (400  ps)  taken  from  refer¬ 

ence  [21].  Note  the  beginning  of  the  collapse  of  the  leading  edge  of  the  pulse  and 
the  series  of  collapse  events  that  consumes  the  pulse  as  it  continues  to  propagate. 

Previous  Work 

As  mentioned  above,  Niday  created  a  code  in  FORTRAN  that  models  the 
propagation  of  a  relatively  long  pulse  (on  the  order  of  picoseconds)  through  air  at 
high  powers,  taking  into  account  nonlinear  terms  that  induce  self-focusing  [21], 
An  example  of  the  collapse  events  and  then  pulse  consumption  resulting  from  the 
propagation  of  a  pulse  in  his  code  is  presented  in  Figure  1.1. 


Niday  came  to  several  conclusions  during  his  investigation  of  the  long  dis¬ 
tance  propagation  of  UV  pulses.  His  first  conclusion  was  that  the  collapses  on 
the  leading  edge  of  the  pulse  are  caused  by  the  transient  shape  of  that  edge.  He 
chose  a  supergaussian  shape  so  that  in  its  middle,  where  it  is  flat,  it  emulates  a 
steady  state  solution  for  which  he  had  solved.  But  the  transient  edge  shape  does 
not  match  the  steady  state  solution,  thus  causing  the  collapses.  He  suggested  that 
the  collapses  are  due  to  the  pulse  trying  to  move  toward  a  steady  state  solution. 
He  also  proposed  that  the  collapses  consumed  the  pulse  at  a  linear  rate,  so  the 
collapses  do  not  develop  faster  as  more  of  the  pulse  is  consumed.  He  was  not, 
however,  able  to  show  convergence  in  time  of  these  collapse  peaks  due  to  limited 
grid  size  and  numerical  parameters.  A  key  component  to  his  research  is  linear 
stability  analysis.  Without  considering  losses,  he  found  that  all  spatial  and  tem¬ 
poral  frequencies  are  able  to  grow,  thus  causing  instability  in  the  propagation  of 
the  pulse. 

Muller  continued  Niday’s  work,  investigating  four  areas  in  further  detail  [8]. 
First,  he  varied  the  numerical  parameters  of  the  propagation  model  in  order  to 
find  the  optimal  combination  of  computational  efficiency  and  resolution  of  pulse 
characteristics  while  striving  for  convergence.  He  concluded  that  a  spatial  grid 
(number  of  data  points  in  the  x  and  y  directions)  of  128  x  128  elements  spaced 
8  fin i  apart  and  a  propagation  step  in  the  z  direction  of  0.5  mm  produced  spatial 
convergence.  In  the  time  domain,  because  of  computational  limitations,  he  sug- 
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gested  a  ratio  of  time  step  to  pulse  duration  of  1:100  for  a  supergaussian  profile 
and  1:20  for  a  Gaussian  on  a  grid  of  512  elements.  This,  however,  did  not  lead  to  a 
fully  converged  result  (see  Figure  1.2).  By  increasing  the  multi-photon  ionization 
(MPI)  coefficient  (3^  by  a  factor  of  10  and  decreasing  the  time  step  to  pulse 
duration  ratio,  he  was  able  to  approach  convergence  for  collapses  of  a  pulse  with 
an  initial  peak  power  of  250  MW  (Figure  1.3).  In  this  work,  time  convergence  of 
the  pulse  collapses  is  produced  using  a  more  powerful  computer  than  what  was 
available  to  Muller. 

Second,  Muller  studied  the  effects  of  initial  power  on  the  collapse  of  the 
pulse.  He  determined  that  higher  initial  power  increased  the  number  of  collapses 
that  occurred,  but  did  not  affect  the  rate  at  which  the  entire  pulse  was  consumed 
by  these  collapse  events  [8].  It  should  be  noted  that  he  did  not  have  complete 
convergence  when  discovering  this  result. 

Third,  Muller  did  a  detailed  study  on  the  relationship  of  the  plasma  to 
the  collapse  events.  He  concluded  that  they  are  very  closely  related.  After  the 
intensity  of  the  pulse  reaches  a  certain  level  due  to  self-focusing,  it  begins  to 
ionize  the  surrounding  air,  thus  generating  plasma.  Once  enough  plasma  is  gener¬ 
ated,  the  pulse  begins  to  defocus,  losing  its  intensity  and  stopping  the  generation 
of  plasma.  As  the  plasma  decays,  self-focusing  begins  to  dominate  once  again, 
creating  another  collapse  event.  He  also  claimed  that  the  nominal  value  of  the 
electron-positive  ion  recombination  coefficient  a  was  1.1  x  10-12  m3/s  since  only 
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(a) 


(b) 


Figure  1.2:  Paul  Muller’s  plot  of  a  500  MW,  20  ps  pulse  after  propagating  25 
cm  for  two  different  values  of  At.  Plot  (b)  represents  a  zoomed  in  area  of  (a), 
and  the  scaling  of  the  plot  has  cut  off  the  top  of  the  first  peak. 
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(a) 


(b) 

Figure  1.3:  Paul  Muller’s  plot  of  a  250  MW,  20  ps  pulse  with  /3^K^  increased  by 
a  factor  of  10.  Plot  (b)  represents  a  zoomed  in  area  of  (a). 
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slight  adjustments  caused  the  model  to  fail.  This  was  discovered  to  be  incorrect 
by  Chalus,  however,  who  suggested  a  much  smaller  value  [24].  This  new  value 
of  1.3  x  10  14  m3/s  is  incorporated  into  this  study.  Muller  also  suggested  that 
collapse  events  produced  on  the  trailing  edge  of  a  pulse  are  due  to  an  overflow  of 
the  model’s  grid  in  the  spatial  frequency  domain.  Care  is  given  in  this  study  to 
alleviate  this  numerical  error. 

Finally,  Muller  added  a  group  velocity  dispersion  (GVD)  term  to  the  propa¬ 
gation  code  to  investigate  much  smaller  pulse  durations.  He  determined  that  the 
GVD  term  did  not  affect  pulses  of  duration  longer  than  100  fs  and  that  even  with 
GVD  taken  into  consideration,  10  fs  pulses  also  collapsed  without  propagating  a 
long  distance.  This  work  will  reproduce  code  to  account  for  GVD  and  incorporate 
it  in  the  study  of  the  collapse  events  even  for  pulses  longer  than  100  fs  because  his 
results  were  obtained  without  full  convergence  of  the  collapse  events.  Considering 
the  collapse  events  are  on  the  order  of  femtoseconds,  GVD  might  play  a  role  in 
their  development. 

Objectives  and  Methods 

This  thesis  continues  the  work  of  Niday,  focusing  specifically  on  Muller’s 
efforts  to  achieve  time  convergence  of  the  collapse  events.  By  updating  the  FOR¬ 
TRAN  code  authored  by  Niday,  the  pulse  duration  is  decreased  from  Niday ’s  400 
ps  pulses  down  to  10  ps  pulses  to  allow  for  a  decrease  in  the  time  grid  spacing 
with  a  reasonable  increase  in  the  number  of  grid  points.  This  provides  enough 
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resolution  to  validate  the  shape  of  the  collapse  events  up  to  a  propagation  distance 
of  ~  50  cm.  Analysis  of  the  time  convergence  and  shape  of  these  collapses  is  done 
with  the  help  of  Matlab®  code  that  Muller  created.  In  order  to  understand  more 
about  what  is  causing  these  collapse  events  and  how  they  are  forming,  a  detailed 
analysis  of  the  individual  terms  in  the  propagation  equation  is  presented  with  the 
help  of  Mathematica®.  All  of  this  is  done  on  a  dedicated  workstation  with  two 
dual-core  processors,  AMD  Opteron  275  with  16  GB  RAM  that  was  procured 
prior  to  the  start  of  this  project. 

An  investigation  of  the  individual  terms  in  the  propagation  equation  is  done 
to  include  the  MPI  coefficient  and  its  role  in  shaping  the  collapse  events.  In 
Muller’s  work,  increasing  this  MPI  parameter  artificially  led  to  better  convergence 
of  the  collapse  events.  An  explanation  for  this  result  is  given  and  discussed.  Also, 
an  updated  value  of  the  electron-positive  ion  recombination  coefficient  a  in  the 
plasma  calculations  is  employed.  Chalus  suggests  that  previous  literature  has 
been  using  a  value  two  orders  of  magnitude  higher  than  the  correct  value  of 
1.3xl0”14  m3/s  [24],  The  difference  this  new  value  makes  on  the  behavior  of  the 
pulse  propagation  is  presented  along  with  an  explanation.  Also,  the  role  that  GVD 
plays  in  the  development  of  collapse  events  is  detailed  to  include  its  small  value 
relative  to  the  other  propagation  terms  and  its  effect  on  the  stability  analysis. 

An  expanded  stability  analysis  based  on  the  analysis  done  by  Niday  in  refer¬ 
ence  [21,23]  is  conducted.  The  analysis  includes  group  velocity  dispersion  effects 
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and  takes  into  account  how  losses  compare  to  the  growth  rate.  The  analysis  is 
done  with  various  combinations  of  the  different  values  for  the  recombination  rate 
and  multi-photon  ionization  coefficient.  It  is  determined  that  the  inclusion  of 
GVD  to  the  diffraction  and  Kerr  terms  that  Niday  used  leads  to  a  slower  growth 
rate  for  high  temporal  frequencies,  except  at  zero  spatial  frequencies  GVD  causes 
the  growth  rate  to  be  non-zero.  Although  a  full  loss  analysis  is  not  conducted, 
the  losses  are  found  to  be  insignificant  compared  to  the  growth  rate,  except  at 
high  intensities  on  the  order  of  101'  W/m2. 


15 


II.  Theory  Review 


This  chapter  will  give  the  reader  a  foundation  in  the  physics  needed  to  understand 
this  research.  It  starts  with  the  basic  principles  that  govern  the  propagation  of 
light,  goes  on  to  address  the  individual  terms  of  the  propagation  equation,  and 
then  goes  into  further  detail  on  the  most  important  terms  for  the  interests  of  this 
thesis. 


Propagation  Equation 

Understanding  the  propagation  of  a  laser  beam  or  pulse  through  air  begins 
with  the  Maxwell  wave  equation, 


n2  d2E  _  1  82PNL 
c2  dt 2  e0c2  dt2 


(2.1) 


where  n0  is  the  linear  refractive  index,  eo  is  the  permittivity  of  free  space,  /x0  is  the 
permeability  of  free  space,  and  c2  =  l/(eo/io)-  Pnl  is  the  nonlinear  polarization 
and  can  include  all  nonlinear  effects,  but  may  be  written  simply  as  6qX^\£\2E 
for  self-focusing,  yU)  is  equal  to  2rion2,  where  n?  is  the  self-focusing  index.  £  is 
a  complex  amplitude  of  the  slowly  varying  envelope  of  the  electric  held  and  has 
units  that  are  defined  such  that  the  irradiance,  or  what  we  will  call  intensity,  \£\2 , 
has  units  of  W/m2.  For  a  wave  propagating  along  the  z-axis,  the  electric  held  E 
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is  represented  by: 


E  (x,  y,  z,  t)  =  if  (x,  y,  z,  t)  e*(t"  “,)  +  c.c. 


(2.2) 


with  uj  being  the  angular  frequency  and  with  k  =  nouj/c.  Using  equations  (2.1) 


and  (2.2)  and  the  nonlinear  polarization  terms  and  effects  of  plasma  generation 


that  Niday  [21,23],  Schwarz  and  Diels  [11],  and  Muller  [8]  use,  one  can  obtain 


to  describe  how  S  changes  as  a  function  of  the  propagation  distance  z.  The  first 
term  on  the  right  hand  side  of  the  equation  represents  transverse  diffraction.  It  is 
followed  by  a  term  expressing  group  velocity  dispersion  (GVD)  where  k"  represents 
the  second  order  derivative  of  the  wave  number  k  with  respect  to  u.  The  third 


term  describes  multi-photon  ionization  where  (5^  is  the  coefficient  for  MPI  of  the 


Kth  order.  In  the  case  of  UV  propagation  in  air,  K  =  3  for  A  =  248  nm  because 
three  photons  with  energy  of  approximately  5  eV  are  required  to  ionize  one  O2 
molecule  with  ionization  energy  of  approximately  13.6  eV,  therefore  generating  a 
plasma.  The  fourth  term  is  a  plasma  term  where  a  is  the  cross  section  for  inverse 
bremsstrahlung,  uj  is  the  optical  reference  frequency  of  the  light,  r  is  the  electron 
collision  time,  and  p  is  the  electron  plasma  density.  The  last  two  terms  represent 
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self-focusing  and  the  Raman  effect,  where  /  is  the  fraction  of  the  Kerr  effect  which 
is  not  instantaneous,  but  delayed  due  to  stimulated  molecular  Raman  scattering. 

Simplifying  Assumptions 

To  be  able  to  model  the  behavior  of  a  propagating  pulse,  taking  into  ac¬ 
count  the  self-focusing  terms  needed  to  produce  long  distance  propagation,  one 
can  make  certain  simplifying  assumptions  in  order  to  accommodate  the  time  and 
computational  abilities  at  hand.  In  equation  (2.3),  we  neglect  the  last  term  de¬ 
scribing  the  Raman  effect  according  to  Niday  [21].  His  investigation  showed  that 
for  UV  pulses  on  the  order  of  picoseconds  the  Raman  contribution  to  n2  was  found 
to  be  insignificant  due  to  the  lack  of  spectral  broadening  before  the  pulses  ionize. 
Niday’s  conclusions  on  Raman  scattering  were  founded  on  the  work  of  Penano  et 
al  [25].  Therefore,  equation  (2.3)  reduces  to 

^  =  4vi  -  ff—\er-2S  -  -  (1  +  *rr)  p£  +  ik0n2\£\>£.  (2.4) 

dz  2k  v±  2  dt2  2  1  1  2  v  011  v  ; 

In  order  to  evaluate  the  plasma  term  in  the  above  equation,  we  describe  the 
electron  plasma  density  p  by 

%  =  Cp|f|2  +  -  Khl  -aP2  +  D  vi  ft  <2'5) 
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where  C  is  the  avalanche  ionization  coefficient,  a  is  the  electron-positive  ion  re¬ 
combination  coefficient,  and  D  is  the  diffusion  strength.  In  previous  literature 
a  value  of  l.lxlO”12  m3/s  was  used  for  a  [5,8,13,21,23].  This  value  was  also 
used  for  the  first  half  of  this  research,  but  was  changed  to  1.3xl0_li  m3/s  to 
agree  with  the  new  value  that  Chalus  presents  in  reference  [24],  We  assume  no 
diffusion,  so  D  =  0,  and  the  pulse  durations  are  short  enough  to  neglect  the 
avalanche  ionization  so  we  can  also  set  C  =  0  [11].  Niday,  Muller,  and  Schwarz 
and  Diels  [8,11,21,23]  use  a  steady  state  solution  found  by  setting  ^  =  0  to 
express  the  plasma  density  as  a  power  of  intensity  for  their  two  dimensional  code. 
We,  however,  only  employ  Niday ’s  3D  code  which  must  solve  the  differential  equa¬ 
tion  (2.5)  for  each  propagation  step.  A  detailed  analysis  of  the  range  needed  for 
these  approximations  is  provided  by  Zhang  [17],  which  says  that  for  steady  state 
to  apply,  the  UV  pulses  must  be  at  least  a  few  tens  of  picoseconds  long,  and  in 
order  to  neglect  avalanche  ionization  they  must  be  shorter  than  several  nanosec¬ 
onds.  Since  we  are  dealing  with  ps  pulses,  the  equation  that  must  be  solved  in 
order  to  describe  the  plasma  density  becomes 


dp  _  /?W \£\2K 
dt  Kfrw 


(2.6) 


For  a  248  nm  pulse  propagating  in  air,  K  =  3,  and  the  value  for  (3^  is  3.9  x  10-34 
m3/W2.  This  value  is  used  in  many  other  references  as  well  [5,8,11,13,21,23]. 
However,  different  values  are  also  employed  in  the  literature  [17,22]  ranging  from 
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9.55  x  10  34  to  1.0  x  10  31 .  This  throws  the  exact  value  of  the  MPI  coefficient  into 


question.  Therefore,  we  will  investigate  what  increasing  this  value  by  an  order  of 
magnitude  does  to  the  propagation  of  the  pulse. 

Another  effect  not  present  in  equation  (2.3)  is  Rayleigh  scattering.  Accord¬ 
ing  to  reference  [17],  after  approximately  400-500  meters  of  propagation,  Rayleigh 
scattering  begins  to  significantly  affect  the  filament  propagation.  The  model  used 
by  Niday  and  Muller  has  not  yet  been  able  to  stably  propagate  to  a  distance  of 
this  magnitude,  but  it  is  worth  mentioning. 


Self-Focusing  and  Filamentation 

A  major  factor  in  filamentation  is  the  self-focusing  Kerr  effect.  The  refrac¬ 
tive  index  is  not  only  frequency  dependent  as  is  addressed  when  dealing  with 
GVD  in  the  following  section,  but  it  is  also  intensity  dependent  according  to 
n  =  no  +  n2l(x,y,t)  [12].  As  the  intensity  increases,  the  nonlinear  Kerr  index 
has  a  greater  effect  on  the  refractive  index.  In  a  normal  Gaussian  pulse,  the  inten¬ 
sity  is  greatest  in  the  center  and  declines  radially  outward.  This,  in  effect,  creates 
a  virtual  lens  that  focuses  the  pulse.  This  self-focusing  can  overcome  diffraction 
once  the  peak  power  of  the  pulse  reaches  a  critical  power  given  by 


P 

J  ri 


A2 

27m0n2  ’ 


(2.7) 
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where  A  is  the  wavelength  and  uq  is  the  refractive  index  of  the  background  [21]. 
This  critical  power  is  125.5  MW  for  a  wavelength  of  248  nm  in  air.  If  allowed  to 
act  in  the  absence  of  other  effects,  self-focusing  can  lead  to  a  catastrophic  collapse 
of  the  pulse,  ffowever,  MPI  and  associated  plasma  defocusing  counteract  this 
effect  to  create  filaments. 

There  are  two  main  models  that  describe  the  phenomenon  of  hlamentation. 
Extensions  of  these  two  models  attempt  to  explain  certain  aspects  of  hlamentation 
with  greater  accuracy.  First,  the  moving  focus  model  considers  the  filament  as 
an  illusion.  The  pulse  is  broken  into  time  slices  that  are  assumed  only  weakly 
coupled  to  each  other.  Each  time  slice  has  its  own  power  which  allows  the  slices 
with  power  greater  than  Pcr  to  be  focused  and  those  with  power  less  than  Pcr  to 
diffract  [12]. 

The  second  model,  self-guiding,  has  two  iterations.  It  begins  with  the  self¬ 
trapping  model  that  suggests  that  a  strict  equilibrium  between  diffraction  and 
self-focusing,  often  called  the  Townes  mode,  is  what  enables  hlamentation  to  oc¬ 
cur.  This  equilibrium  was  adjusted  with  the  self-channeling  model  to  loosen  the 
strictness  of  the  balance  and  include  refraction  from  plasma  generated  by  the  in¬ 
tense  laser  pulse.  This  addressed  the  problem  of  instability  of  self-trapping  that 
led  to  either  an  inevitable  collapse  of  high  powered  pulses  or  an  eventual  diffrac¬ 
tion  of  low  powered  pulses  on  the  same  order  as  normally  propagated  beams  [12]. 
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The  dynamic  spatial  replenishment  model  is  a  combination  of  the  previous 
two  main  models.  This  model  shows  that  the  leading  edge  of  the  pulse  forms  a 
peak  due  to  self-focusing  which  generates  a  plasma  that  causes  the  trailing  edge 
of  the  pulse  to  defocus  and  the  leading  peak  to  lessen  in  intensity  due  to  MPI.  As 
the  peak  intensity  is  lowered,  the  plasma  is  then  turned  off,  allowing  the  pulse  to 
self-focus  once  again  until  this  scenario  repeats  itself  enough  to  make  the  pulse 
power  insufficient  for  refocusing  [12].  This  model  allows  for  the  long  distance 
propagation  of  the  filaments  and  is  the  foundation  for  this  research. 

The  spatiotemporal  soliton  or  light  bullet  is  the  proposed  idea  of  a  soliton 
forming  in  both  space  and  time  due  to  nonlinear  spatial  effects,  such  as  self- 
focusing,  that  compensate  for  diffraction,  and  self-phase  modulation  that  com¬ 
pensates  for  the  spreading  of  the  beam  in  time.  However,  this  proposal  is  highly 
unlikely  and  has  never  been  observed  experimentally.  This  is  probably  due  to  the 
instability  of  the  Townes  modes  discussed  earlier  and  the  neglect  of  physical  ef¬ 
fects  always  present  in  practice  that  would  prevent  a  light  bullet  from  propagating 
longer  distances  than  normal  diffraction  and  dispersion  lengths  [12]. 

X-waves  or  Bessel-X  waves  are  one  more  suggested  solution  to  explain  fil¬ 
ament  at  ion .  They  are  wave  packets  that  propagate  without  undergoing  natural 
spreading  due  to  diffraction  or  dispersion.  They  form  an  X  shape  in  both  the 
near  and  far  fields  as  represented  in  Figure  2.1.  These  peculiar  waves  come  from 
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a  localized  stationary  solution  of  the  form 


£(r,  t) 


1 

2tt 


du)  f  (lu)  Jq{V  kk" 


\u>\r)exp(—icut) , 


(2.8) 


for  equation  (2.4)  if  only  the  first  two  right  hand  side  terms  are  considered.  In 
equation  (2.8)  /( u)  is  any  narrow  spectral  function  [12],  Because  X-waves  have 
been  observed  to  share  many  characteristics  of  femtosecond  filaments,  they  are  a 


promising  concept  for  the  interpretation  of  hlamentation. 
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Figure  2.1:  This  figure  shows  how  an  X-wave  has  the  same  shape  in  both  (a) 
the  r-t  domain  or  near  field  and  (b)  the  k-u)  domain  or  far  field.  Taken  from 
Ref  [12], 


Group  Velocity  Dispersion 

There  are  many  different  effects  that  can  be  included  in  an  investigation  of 
the  long  distance  propagation  of  laser  pulses,  depending  on  the  assumptions  and 
approximations  made.  Niday  [21,23]  and  Schwarz  and  Diels  [11]  did  not  include 
the  GVD  term  in  their  propagation  equation  due  to  the  assumption  that  the  pulses 
they  were  dealing  with  were  long  enough  to  ignore  GVD.  This  assumption  was 
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reasonable  for  the  initial  pulse  durations,  but  with  the  collapse  events  generated 


in  Niday’s  code  being  on  the  order  of  femto-  and  pico-  seconds,  it  is  important 


to  revisit  the  relevance  of  GVD.  Muller  claimed  that  GVD  did  not  affect  pulses 


longer  than  100  fs  in  duration  [8],  but  since  he  could  not  achieve  convergence  of 
the  collapses,  he  was  not  able  to  confirm  that  GVD  did  not  effect  the  formation  of 


these  collapses.  Although  the  pulses  being  studied  here  are  long  (on  the  order  of 


picoseconds),  the  group  velocity  dispersion  term  must  not  be  ignored  due  to  the 
collapse  events  having  durations  on  the  order  of  tens  of  femtoseconds.  With  parts 
of  the  propagating  pulse  being  so  short,  GVD  should  be  included  to  accurately 
model  the  collapse  event’s  behavior  during  propagation. 

GVD  comes  from  the  dependence  that  a  wave’s  velocity  has  on  its  frequency. 
This  causes  the  different  frequencies  of  the  wave  to  spread  out  temporally  as  a 
result  of  the  different  velocities.  It  works  in  a  similar  way  to  how  diffraction 
spreads  a  pulse  out  spatially.  Extracting  the  GVD  portion  from  our  propagation 
equation  (2.3)  we  have 


d£_  _  ik"  d2£ 
dz  2  dt2 ' 


(2.9) 


With  a  little  rearranging,  this  can  then  be  compared  to  the  one  dimensional  spatial 
propagation  of  a  Gaussian  beam  as  described  by 


(2.10) 
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from  equation  (4.7.11)  in  reference  [26].  Thus,  to  solve  for  the  GVD  parameter, 
one  can  use  a  similar  process  as  with  a  ID  Gaussian  beam. 

In  order  to  ensure  accuracy,  the  value  of  k"  was  explicitly  calculated,  since 
there  are  varying  values  in  the  literature.  Recognizing  the  difference  in  field 
representations  of  equations  (2.2)  and  Verdeyen’s  (4.7.4)  with  opposite  signs  in 
the  exponentials,  one  concludes  that  Verdeyen’s  /?2  is  equivalent  to  k"  in  equation 
(2.9).  Verdeyen  goes  on  to  define  02  as  (A3/27tc2)  •  (d2n/d\2).  To  determine  this 
value,  one  needs  the  index  of  refraction  for  air  in  terms  of  wavelength.  This  is 
obtained  using  Cauchy’s  Formula  [27] 


n„ir  =  1.000287566  + 


1.34  x  10_18m2  3.777  x  10' 

—  +  — 


m 


A2 


A4 


(2.11) 


Using  this  information,  a  value  of  k"  =  1.45  fs2/cm  was  obtained  for  A  =  248 
nm.  This  value  is  what  is  used  in  this  study  and  is  between  Couairon’s  value 
of  k"  =1.21  fs2/cm  [9]  and  Schwarz’s  value  of  k"  =2  fs2/cm  [22],  The  fact  that 
this  term  is  positive  affirms  that  we  are  dealing  with  normal  and  not  anomalous 
dispersion. 


The  reason  GVD  plays  a  part  in  the  broadening  of  the  duration  of  these 
collapse  events  is  the  dependence  of  the  dispersion  length,  Lp,  on  the  pulse  half¬ 
width  in  time,  tp,  of 


Ld  — 


t2 

Lp 


k"' 


(2.12) 
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If  the  pulse  duration  is  small,  this  dispersion  length,  or  distance  the  pulse  travels 
before  it  is  temporally  broadened  by  a  factor  of  y/2,  will  also  be  small  allowing 
GVD  to  broaden  the  pulse  in  the  distances  being  propagated  during  these  simu¬ 
lations.  In  accordance  with  the  widths  obtained  for  the  collapse  events  once  the 
collapses  form,  Lp  ra  20  cm.  By  including  this  GVD  parameter  in  the  code,  it  is 
thought  that  the  collapsing  events  of  the  pulse  may  be  limited  by  the  temporal 
dispersion  effects. 

Theoretical  Summary 

In  summary,  the  propagation  equation  used  to  describe  the  behavior  of  a 
pulse  of  light  is  based  on  Maxwell’s  equations.  The  terms  included  in  this  study 
represent  transverse  diffraction,  multi-photon  ionization,  plasma  generation  and 
absorption,  Kerr  self-focusing,  and  group  velocity  dispersion.  This  last  term  was 
added  to  the  study  because  of  the  small  duration  of  the  collapse  events  being 
investigated.  The  value  of  k"  for  GVD  was  established  as  1.45  fs2/cm.  Terms 
such  as  the  Raman  effect,  Raleigh  scattering,  avalanche  ionization,  and  plasma 
diffusion  were  neglected  according  to  valid  assumptions  based  on  the  research 
parameters.  Long  distance  propagation  is  obtained  by  hlamentation  of  the  pulse 
when  the  power  reaches  a  certain  critical  power  to  induce  self-focusing.  Many 
theoretical  models  have  been  proposed  to  explain  this  hlamentation  such  as  the 
dynamic  spatial  replenishment  model  which  suggests  somewhat  of  a  tug-of-war 
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game  on  the  intensity  by  self-focusing  and  MPI  that  keeps  the  pulse  propagating. 
There  is  also  the  possibility  that  X-waves  could  be  a  description  of  hlamentation. 
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III.  Methodology 


This  chapter  will  describe  the  basic  physical  and  numerical  parameters  used  in  the 
simulations  of  the  propagation  of  UV  pulses.  A  table  presenting  the  parameters 
used  in  the  bulk  of  the  simulated  runs  is  provided  for  the  convenience  of  the  reader. 
We  will  go  into  the  details  of  how  the  code  is  written,  to  include  an  explanation 
of  the  split  operator  method  and  what  limitations  are  incurred  due  to  numerical 
analysis.  Details  of  the  modification  of  the  original  code  of  Niday’s  will  also  be 
discussed. 

Physical  and  Numerical  Parameters 

A  summary  of  the  physical  and  numerical  parameters  used  in  this  study 
will  be  very  useful  to  the  reader  when  trying  to  understand  exactly  how  the 
results  presented  in  this  thesis  were  obtained.  The  key  input  parameters  for  the 
initial  condition  of  the  pulse  are  waist  size,  pulse  duration,  pulse  power  or  peak 
intensity,  and  desired  grid  size.  The  default  initial  parameters  used  are  specified 
in  Table  3.1.  A  few  of  the  actual  grid  parameters  varied  quite  a  bit  from  run 
to  run  such  as  Nt  and  At,  so  they  are  not  included  in  this  table.  They  will  be 
specified  appropriately  when  data  is  presented  along  with  any  parameters  that 
are  different  from  the  values  stated  here. 

A  few  notes  on  these  values  are  needed  to  ensure  clarity.  The  values  of 
Nx,  Ny,  Ax  and  Ay  were  used  based  on  the  parameter  study  by  Muller  [8]  for 
the  waist  size  Wq.  The  supergaussian  order,  M,  is  used  in  the  time  dimension  to 


Table  3.1:  Parameters  for  UV  pulse  simulations 


Parameter 

Variable 

Value 

Units 

Number  of  grid  points  in  x-direction 

Nx 

128 

— 

Spacing  of  grid  points  in  x-direction 

Ax 

8 

/irn 

Number  of  grid  points  in  y-direction 

Ny 

128 

— 

Spacing  of  grid  points  in  y-direction 

Ay 

8 

/irn 

Wavelength 

A 

248 

nm 

Supergaussian  order 

M 

10 

— 

Pulse  Duration  (half-width) 

tp 

5 

ps 

Distance  of  each  propagation  step 

A  z 

0.5 

mm 

Input  power 

Pin 

250 

MW 

Critical  power 

P 

1  cnt 

125.5 

MW 

Pulse  waist 

w  o 

120 

/irn 

Nonlinear  index  of  refraction 

n2 

7.8xl0-23 

m2/W 

Linear  loss 

OIL 

2.5x  10~4 

m”1 

MPI  order 

K 

3 

— 

MPI  coefficient 

3.9xl0-34 

m3/W2 

Inverse  Bremsstrahlung  cross  section 

a 

5.2xl0”25 

m2 

Electron  collision  time 

T 

0.35xl0~13 

ps 

Old  radiative  recombination  rate 

a 

l.lxlO-12 

m3/s 

New  radiative  recombination  rate 

1.3xl0-14 

m3/s 

GVD  Parameter 

k" 

1.45154xl0-28 

s2/m 
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simulate  the  beginning,  middle,  and  end  of  an  actual  pulse  with  duration  tp  as 
defined  by  a  one  dimensional  Gaussian  of 

( t  \2M 

E(t)  =  E0e  .  (3.1) 

This  is  in  order  to  simulate  a  steady  state  for  much  of  the  duration  of  the  pulse. 
However,  the  beginning  and  the  end  of  the  pulse  will  not  satisfy  the  steady  state 
solution.  The  old  radiative  recombination  rate  in  the  table  was  initially  used  based 
on  many  references  [5,8,13,21,23].  This  value  was  supposedly  an  improvement  on  a 
value  used  by  Schwarz  of  5x  10”13  m3/s  [22],  Recently,  though,  Chalus  presented  a 
newer  value  for  a  which  is  defined  in  the  table  as  the  ‘new’  radiative  recombination 
rate  an.  Many  of  the  simulations  done  in  this  research  were  performed  with  this 
updated  value  and  will  be  specified  when  appropriate. 

Numerical  Methods 

All  simulations  done  for  this  research  used  FORTRAN  code  derived  from 
what  was  originally  written  by  Niday  [21],  His  three  dimensional  code  uses  a  three 
dimensional  grid  or  array  to  hold  the  information  representing  the  field.  Symmetry 
in  the  x  and  y  directions  is  assumed  for  the  initial  condition  of  the  Gaussian  pulse. 
Time  is  the  third  dimension  of  the  array  and  can  be  related  to  the  third  spatial 
dimension  by  z  =  cAt  where  c  is  the  speed  of  light,  At  is  the  amount  of  time  that 
has  passed,  and  z  is  the  direction  of  beam  propagation.  At  each  propagation  or 
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step,  the  pulse  is  represented  with  an  x,  y,  and  t  array  where  t  —  0  is  located 
at  z.  The  reference  frame  moves  at  a  speed  equivalent  to  the  group  velocity  of  the 
pulse  in  order  to  record  the  pulse’s  evolution  at  each  distance  z.  If  the  reference 
frame  was  fixed,  the  pulse  would  go  off  the  time  grid  almost  immediately  after 
beginning  to  propagate  because  the  speed  of  light  is  so  much  larger  than  the  size 
scale  of  the  pulse  and  grid.  By  allowing  the  reference  frame  to  move  with  the 
pulse,  we  are  able  to  use  grid  sizes  and  spacings  that  can  accurately  describe  the 
pulse  as  it  propagates.  As  the  grid  spacing  decreases,  the  grid  size  (number  of 
points  on  the  grid)  must  increase  to  account  for  the  same  sized  pulse.  This,  in 
turn,  dramatically  increases  the  amount  of  time  and  memory  space  needed  to 
compute  the  information.  Thus,  there  is  a  delicate  balance  between  resolution 
and  practicality.  Figures  3.1  through  3.3  are  an  example  of  the  3D  plots  that  can 
be  generated  by  the  GUI  interface  associated  with  Niday’s  code.  These  figures 
show  the  evolution  of  a  pulse  with  a  time  history  plot  of  the  field  with  respect 
to  the  x  axis  along  with  its  fluence  profile  at  various  propagation  distances.  It  is 
a  ten  picosecond  pulse  generated  with  the  new  recombination  rate  and  no  GVD, 
represented  on  a  grid  with  Nt  =  4096  and  At  =  3  fs.  Note  that  at  50  cm  the 
fluence  profile  begins  to  look  square.  This  means  it  is  no  longer  an  accurate 
representation  since  the  fluence  has  gone  off  the  spatial  grid  and  has  begun  to  be 
confined  by  the  grid,  the  absorbing  boundary  conditions  reflecting  nonzero  values. 
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In  Niday’s  3D  propagation  code,  the  x,  y,  and  t  dimensions  are  transformed 
via  a  Fast  Fourier  Transform  (FFT)  as  required  by  the  split,  operator  method 
which  is  explained  in  the  following  paragraph.  In  order  for  the  FFT  to  work 
efficiently,  the  number  of  grid  points  in  each  dimension  need  to  be  chosen  as 
powers  of  two.  Originally,  Niday’s  code  contained  two  3D  arrays  to  represent  the 
held:  one  array  for  the  x,  y,  and  t  dimensions  in  real  space  and  the  other  array  for 
the  kx,  ky,  and  D  dimensions  in  frequency  or  k-  space  after  the  Fourier  transform. 
Having  a  copy  of  the  held  in  both  the  space/time  and  frequency  domains  allowed 
for  more  freedom  in  the  manipulation  of  the  held  within  the  code.  However,  as 
resolution  needed  to  be  increased,  the  size  of  the  grid  also  had  to  be  increased. 
With  the  original  code,  when  the  grid  size  was  doubled,  the  memory  needed  to 
hold  both  of  the  held  arrays  quadrupled.  By  eliminating  the  extra  array  dedicated 
to  the  fc-space  representation  of  the  held,  the  grid  size  was  able  to  double  once 
more  beyond  what  was  previously  allowed.  To  do  this,  the  code  was  adjusted 
to  keep  only  one  array  in  memory  for  the  held.  This  array,  therefore,  switched 
between  real  and  fc-space.  With  the  ability  to  have  a  larger  grid  size,  the  resolution 
in  the  time  domain  was  able  to  be  increased.  This  helped  to  converge  the  collapse 
events  in  the  time  domain. 

The  split  operator  method  is  used  by  Niday  in  his  code  to  solve  the  pulse 
propagation  equation  (2.4),  which  can  be  simplified  to 

^  =  iTS  +  iVS , 
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(3.2) 


z=12.5  cm 


^  (picoseconds) 
(c) 


Figure  3.1:  Left  Column:  Full  3D  view  of  pulse  propagation  with  its  leading 

edge  to  the  right  in  the  negative  time  domain.  Right  Column:  Corresponding 
fluence  profiles  at  propagation  distances  of  0  and  12.5  cm. 
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Figure  3.2:  Left  Column:  Full  3D  view  of  pulse  propagation  with  its  leading 

edge  to  the  right  in  the  negative  time  domain.  Right  Column:  Corresponding 
fluence  prohles  at  propagation  distances  of  25  cm  and  50  cm.  Note  the  evidence 
in  the  Auence  proAle  of  going  oA  the  spatial  grid. 
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z=75  cm 


Figure  3.3:  Left  Column:  Full  3D  view  of  pulse  propagation  with  its  leading 

edge  to  the  right  in  the  negative  time  domain.  Right  Column:  Corresponding 
fluence  profiles  at  propagation  distances  of  75  cm  and  1  m.  Figure  3.4  shows  that 
the  pulse  has  not  gone  off  the  k  and  grids  by  these  propagation  distances. 
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where  we  have 


iV  =  ik0n2\S\2  — 


2  P^K\c\2K-2  a 


-w 


—  —  (1  +  iur)  p  and 


(3.3) 


-  .  1  o  ik"  d 2 

'  =  ‘2ki, v±  2 "ae  ' 


(3.4) 


Equation  (3.3)  is  the  nonlinear  index  change  piece  and  will  be  diagonal  in  the 
spatial  representation  for  £.  Equation  (3.4)  is  the  linear  diffraction  and  GVD 
piece  and  is  diagonal  in  the  spectral  representation  of  £  or  in  fc-space.  The  split 
operator  method  represents  the  incremental  step  propagator  elHAz 


as 


•  1  LAZ  •  T  7"  A  •  1  LAZ  ,  0s 

e.—  e,VAze,—  +  0(A£3), 


(3.5) 


where  iH  =  if  +  iV.  In  order  to  use  this  method  effectively  to  solve  for  the 
solution  £(z  +  Az)  =  elHAz£(z ),  we  apply  the  three  separate  exponential  parts 
of  equation  (3.5)  at  different  times  while  transforming  £  between  its  real  and 
spectral  representations.  First,  the  FFT  is  taken,  putting  £  into  fc-space,  then 
the  first  half  of  the  T  operator  is  applied  in  its  spectral  representation  according 
to 


'  1  LAZ  /  - 

e  2  S(kx,  ky)  =  e 


k 

4fc0 


i  k  n- 


Az 


ky 


(3.6) 


This  spectral  representation  of  the  T  operator  can  be  determined  for  each  point 
on  the  k  grid,  that  is,  each  value  of  k±  and  at  the  beginning  of  the  code 
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after  the  first  FFT.  It  does  not  have  to  be  recomputed  for  each  propagation  step. 


Transforming  equation  (3.6)  back  into  its  spatial  representation  with  an  inverse 
FFT,  the  nonlinear  index  change  piece  is  then  applied  by  multiplying  by  elVAz. 
Another  Fourier  transform  puts  £  back  into  spectral  space  to  allow  for  the  second 
half  of  the  T  operator  to  be  applied  in  the  same  way  as  the  first  half. 

Limitations  of  Numerical  Methods 

Using  this  split  operator  method  creates  a  requirement  for  the  grid  size  and 
spacing  to  ensure  accuracy.  Both  the  spatial  and  spectral  grids  must  contain 
adequate  resolution  and  extent  to  represent  the  full  held.  For  this  split  operator 
approximation  to  work  well,  the  values  in  the  exponentials  should  be  small  so 
that  the  arguments  of  the  exponential  phase  factors  will  be  smaller  than  unity. 
Including  the  GVD  term,  we  have  for  the  maximum  absolute  value  of  T, 

k\  k"£l2 

max  1 7)  I  max 

Zh^Q  Z 

The  value  on  the  right  hand  side  is  maximized  when  the  bigger  of  the  k_\_  and 
terms  is  at  its  max  and  the  smaller  is  zero  since  both  terms  are  always  positive. 
For  Ax  =8  /irn,  At  has  to  be  less  than  0.5  fs  to  allow  the  term  to  dominate. 
Since  the  smallest  At  we  use  is  0.9  fs,  we  set  fl2  equal  to  zero  and  have  just  the 


(3.7) 


37 


k±  term  with  which  to  deal.  This  gives  a  maximum  exponential  argument  of 


\f\maxAz  _  n2Az 
2  “  4k0Ax2  ' 


Setting  this  to  be  less  than  unity  and  solving  for  Az  gives 


Ax;  < 


Ak0Ax2 


(3.8) 


(3.9) 


which  is  the  same  condition  that  Niday  found  [21].  With  an  UV  wavelength  of 
248  nm  and  Ax  =8  fin i,  Ax;  <  0.67  mm.  This  condition  is  satisfied  with  our 
choice  of  Ax;  =  0.5  mm  (see  Table  3.1).  If  At  is  small  enough  to  allow  the  Q  term 
to  dominate  in  equation  (3.7),  Ax;  would  then  be  limited  by 


Ax;  < 


4At2 

k"  it2 


(3.10) 


The  biggest  limitations  of  numerically  solving  the  propagation  equations 
to  describe  UV  hlamentation  are  computing  power  and  time.  A  computer  with 
two  dual  processors  totaling  16  GB  of  RAM  has  been  used  for  the  computations 
presented  here.  The  longest,  most  detailed  runs  took  up  to  six  days  to  complete. 
The  smaller  runs  (propagating  less  distance  and  with  less  resolution)  took  12 
to  36  hours,  thus  limiting  how  many  runs  were  able  to  be  performed  and  the 
resolution  able  to  be  achieved  in  the  time  frame  given  for  this  research.  The  code 
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was  modified  to  use  almost  50%  less  memory  in  order  to  be  able  to  achieve  a  time 
step  of  0.9  fs  that  was  needed  to  show  convergence.  Since  this  run  took  over  12 
GB  of  the  16  total  available  and  the  FFT  requires  the  dimensions  of  the  grid  to  be 
in  powers  of  2,  increasing  the  grid  size  to  accommodate  a  smaller  change  in  time 
would  require  almost  a  doubling  of  the  memory  needed  for  computation,  putting 
it  well  above  the  computer’s  limit. 

Another  facet  of  numerical  limitation  involves  the  fc-space  representation  of 
the  pulse.  In  Muller’s  research,  he  found  that  after  a  certain  distance  of  propa¬ 
gation,  the  trailing  edge  of  the  pulse  would  start  to  exhibit  collapse.  This  was 
thought  to  be  due  to  numerical  parameters  chosen  that  allowed  the  field  to  over¬ 
flow  the  Q-grid,  causing  the  values  to  wrap  around  the  grid.  This  research,  how¬ 
ever,  does  not  see  this  manifestation  of  the  grid  overflow.  Figure  3.4  shows  the 
k  representation  of  the  pulse  in  Figures  3.1  through  3.3  at  propagation  distances 
of  2.5  cm,  12.5  cm,  50  cm,  and  1  m.  Notice  how  the  pulse  does  not  go  off  the  k- 
or  O-grid,  thus,  in  Figures  3.2  and  3.3  the  pulse  does  not  produce  collapses  on  its 
trailing  edge.  Also  notice  how  the  fluence  profile  begins  to  show  the  squareness  of 
the  grid  rather  than  being  spherical  at  only  50  cm.  This  means  that  the  pulse  has 
diverged  to  the  extent  of  going  off  the  spatial  grid.  The  results  of  the  propagation 
code  after  this  point  are  no  longer  accurate  due  to  the  grid’s  artificial  limitations. 


39 


(c)  (d) 

Figure  3.4:  Frequency  representation  of  the  propagating  pulse  from  Figures  3.1 
through  3.3  at  2.5  cm,  12.5  cm,  50  cm,  and  1  m.  Note  that  the  pulse  has  not  gone 
off  the  kx  or  hi  grid.  All  plots  are  on  the  same  scale. 
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To  help  prevent  this  wrap  around  effect  in  the  space  dimensions,  absorbing 


boundary  conditions  (ABC’s)  are  used  in  the  code  according  to 


^Texp{ - w - ) 

SA-  (  (f(x,y)-f(Xmax,y))2\ 

—4—  exP  y - W - ) 

SA-  (  (f(x,y)~  f(x,ymin))2\ 

^Texp[ - w - ) 

SA-  (  (f(x,y)  ~  f(x,ymax))2\ 

^Texp[ - w - )’ 


(3.11) 


where  S  is  the  strength  of  the  ABC’s  of  5000  and  W  is  a  width  of  0.5  x  10-8. 
The  field  is  multiplied  by  eabc  during  the  middle  step  of  the  split  operator  method 
while  in  real  space.  This  subdues  the  values  of  the  field  at  the  edges  of  the  spatial 
grid,  preventing  them  from  being  wrapped  around  to  the  other  side  of  the  grid, 
while  making  sure  not  to  eliminate  values  big  enough  to  make  a  difference  in  the 
pulse  profile.  ABC’s  were  not,  however,  applied  in  the  time  dimension  for  fear 
that  important  information  relevant  to  the  propagation  of  the  pulse  would  be 
ignored.  This  required  the  time  grid  to  maintain  a  large  buffer  between  the  edge 
of  the  pulse  and  the  edge  of  the  grid  to  allow  for  dispersion. 


Summary  of  Methodology 

In  summary,  the  code  written  by  Niday  in  FORTRAN  to  describe  the  prop¬ 
agation  of  a  pulse  through  air  uses  the  split  operator  method  to  apply  a  linear 
operator  to  the  field  in  real  space  and  a  nonlinear  operator  in  k- space  by  using  a 
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FFT  to  convert  the  field.  ABC’s  were  applied  to  prevent  diffracting  field  values  at 
the  edge  of  the  spatial  grid  from  wrapping  around  and  causing  nonphysical  effects 
in  the  propagation  code.  Using  the  split  operator  method  creates  a  requirement 
for  the  grid  parameters,  which,  in  this  case,  restricts  the  z-step  to  an  upper  limit 
of  6.7  mm.  The  numerical  parameters  are  also  limited  by  computational  power 
and  efficiency.  The  code  was  modified  to  use  less  memory  in  order  to  increase  the 
grid  size  to  Nt  =  16384  and  decrease  the  grid  spacing  to  At  =  0.9  fs  in  the  time 
domain  to  achieve  maximum  resolution  with  the  equipment  available.  The  values 
of  most  of  the  parameters  are  specified  in  Table  3.1. 
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IV.  Results  and  Analysis 


This  chapter  includes  the  important  results  obtained  from  the  research  and  simu¬ 
lations  dealing  with  the  long  distance  propagation  of  UV  light.  It  shows  that  time 
convergence  was  achieved  with  certain  parameters  and  explains  why  convergence 
could  not  be  achieved  at  propagation  distances  farther  than  specified.  It  discusses 
the  details  of  the  collapse  events  obtained  from  the  convergence  runs  as  well  as  the 
interaction  between  the  different  terms  of  the  propagation  equation,  explaining  at 
what  points  certain  terms  dominate  and  act  to  increase  or  decrease  this  collapse 
phenomenon.  A  stability  analysis  is  done  with  the  inclusion  of  GVD,  the  new 
recombination  coefficient,  and  the  higher  MPI  parameter.  Thought  was  given  to 
the  effects  of  the  loss  terms  as  well.  In  all  cases,  the  analysis  shows  the  inherent 
instability  of  the  pulses  being  considered.  Suggestions  are  given  to  create  a  more 
stable  propagation.  The  results  of  a  simulation,  using  the  same  initial  conditions 
as  an  experiment  producing  UV  hlamentation,  are  also  presented. 

Time  Convergence 

In  order  to  understand  more  about  these  collapse  events  on  the  leading  edge 
of  the  pulse,  time  convergence  must  be  achieved  with  the  numerical  parameters 
to  make  sure  that  the  information  being  produced  is  accurately  displaying  what 
is  physically  happening.  Using  a  pulse  of  100  ps  total  duration  (tp  =  50  ps),  time 
convergence  was  achieved  with  an  increased  MPI  coefficient,  (5^K)  =  3.9  x  10~33 
m3/W2,  and  no  GVD  with  a  At  of  20  fs.  Figure  4.1  shows  this  convergence  at  a 
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propagation  distance  of  1  m.  The  reason  higher  MPI  leads  to  better  convergence 
of  the  collapses  is  that  it  increases  the  losses  enabling  them  to  better  counteract 
self-focusing.  With  self-focusing  abated  more  quickly,  the  collapses  do  not  have 
an  opportunity  to  increase  to  such  a  great  extent  in  intensity  or  become  so  short 
in  duration,  thus  allowing  easier  representation  with  only  a  sampling  of  the  data. 
With  the  original  value  of  the  MPI  coefficient,  the  intensities  of  the  collapse  events 
are  over  an  order  of  magnitude  larger  than  what  is  shown  here.  This  method  of 
increasing  (3^  was  the  only  way  Muller  could  achieve  convergence.  In  this  effort, 
convergence  was  also  achieved  without  the  artificially  increased  parameter  by 
shortening  the  pulse  to  10  ps,  so  that  the  numerical  parameters  could  be  tightened 
while  the  computer’s  processing  capabilities  could  still  manage  the  data. 

In  Figure  4.2,  a  pulse  consisting  of  the  old  radiative  recombination  rate  of 
1.1  x  1CV12  m3/s,  Nt  =  8192,  and  At  =  1.35  fs,  in  addition  to  the  parameters 
in  Table  3.1  shows  convergence  at  a  propagation  distance  of  20  cm  as  it  is  repre¬ 
sented  in  exactly  the  same  fashion  as  the  same  pulse  with  even  tighter  numerical 
parameters  of  Nt  =  16,  384  and  At=0.9  fs.  However,  going  beyond  this  propaga¬ 
tion  distance,  the  GVD  causes  the  pulse  to  break  down,  which  will  be  discussed  in 
a  later  section.  To  check  the  convergence  at  farther  propagation  distances,  runs 
were  made  with  no  GVD  included.  As  the  propagation  distance  increased,  the 
collapse  events  continue  to  decrease  in  duration  causing  convergence  to  become 
more  difficult.  Figure  4.3  shows  successful  convergence  at  z  =  25  cm  and  slightly 
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Figure  4.1:  This  time  slice  at  x  =  0  and  y  =  0  shows  time  convergence  of  a 

100  ps  pulse  at  z  =  1  m.  The  initial  peak  amplitude  is  normalized  to  1.  Plot  (b) 
is  a  magnification  of  the  collapses  in  plot  (a).  Note  that  there  are  actually  three 
different  lines  on  the  graph,  but  convergence  makes  them  indistinguishable. 
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Figure  4.2:  This  plot  shows  time  convergence  of  the  pulse  collapse  at  z  =  20 
cm.  Plot  (b)  is  a  magnified  portion  of  the  collapse  in  plot  (a).  Note  that  the  two 
lines  are  barely  distinguishable  due  to  convergence. 


46 


less  successful  convergence  at  z  =  50  cm.  It  is  interesting  to  note  that  the  first 
collapse  event  is  always  accurately  displayed  in  duration,  even  if  the  sampling 
frequency  does  not  allow  accurate  display  of  the  peak  intensity  of  the  collapse. 
The  following  collapses  are  shown  to  occur  slightly  later  as  At  is  decreased.  This 
is  believed  to  be  caused  by  the  amount  of  power  in  the  first  collapse  calculated  by 
the  numerical  simulation.  As  At  decreases,  the  volume  under  the  first  collapse  is 
more  accurately  assessed  as  being  larger.  With  the  rest  of  the  pulse  having  less 
power  with  which  to  work,  a  longer  time  is  required  for  a  new  collapse  event  to 
form,  thus  explaining  the  discrepancy  in  the  placement  of  the  following  collapses. 

Time  convergence  was  also  displayed  with  the  new  recombination  rate,  an 
and  no  GVD.  See  Figure  4.4.  A  comparison  was  done  between  the  old  and  new 
recombination  rates  to  see  how  this  value  affected  the  pulse  formation.  As  shown 
in  Figure  4.5,  the  value  of  a  did  not  affect  the  formation  of  the  first  collapse  at 
all.  However,  the  new  smaller  rate  required  more  time  before  the  formation  of 
the  second  collapse  could  occur.  This  is  intuitive  when  one  considers  equation 
(2.6).  A  higher  a  means  slower  change  in  the  plasma  density.  Therefore,  it  takes 
longer  for  the  plasma  to  decrease  enough  for  self-focusing  to  begin  forming  the 
next  collapse.  A  plot  of  the  plasma  densities  comparing  both  a  values  at  z  —  12.5 
cm  and  z  —  25  cm  is  shown  in  Figures  4.6  and  4.7  along  with  the  pulse  profiles 
at  the  same  distances. 
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(a)  z  =  25  cm 


(b)  First  Collapse 


(c)  z  =  50  cm 


(d)  First  Collapse 


(e)  Second  Collapse 


Figure  4.3:  This  plot  shows  time  convergence  of  the  pulse  collapse  without 

GVD  and  a  recombination  rate  of  a  =  1.1  x  10~12  m3/s  at  z  —  0.25  m  (a,b)  and 
z  =  50  cm  (c,d,e).  Notice  that  at  50  cm,  the  peak  of  the  first  collapse  (d)  is  being 
cut  off  due  to  sampling  frequency,  thus  making  the  second  collapse  (e)  slightly 
different  than  the  smaller  time  step  collapse. 
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Figure  4.4:  This  plot  shows  time  convergence  of  the  pulse  collapse  without 

GVD  and  a  recombination  rate  of  an  =  1.3  x  ICC11  m3/s  at  z  —  0.25  m  (a,b)  and 
z  =  50  cm  (c,d,e).  Notice  that  at  50  cm,  the  peak  of  the  first  collapse  (d)  is  being 
cut  off  due  to  sampling  frequency,  thus  making  the  second  collapse  (e)  slightly 
less  converged. 
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Figure  4.5:  A  pulse  with  both  the  old  (solid)  and  new  (dotted)  recombination 
rates  is  shown  at  z  =  50  cm  with  At  =  0.9  fs  and  without  GVD.  (b)  Both  rates 
form  exactly  the  same  first  collapse,  (c)  The  second  collapse  forms  much  later 
with  the  new  recombination  rate. 
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Figure  4.6:  (a)  Illustrates  how  the  plasma  density  with  an  (dotted  line)  de¬ 

creases  at  a  much  slower  rate  than  the  plasma  density  with  a  (solid  line).  Shown 
for  z  =  12.5  cm  and  no  GVD.  (b)  Corresponding  pulse  at  z  =  12.5  cm. 
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Figure  4.7:  (a)  Illustrates  how  the  plasma  density  with  an  (dotted  line)  de¬ 

creases  at  a  much  slower  rate  than  the  plasma  density  with  a  (solid  line).  Shown 
for  z  =  25  cm  and  no  GVD.  (b)  Corresponding  pulse  at  z  =  25  cm. 
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As  the  collapses  continue  to  propagate,  they  become  increasingly  intense  and 
narrow.  This  made  it  difficult  to  assess  the  behavior  of  these  collapses  at  farther 
propagation  distances.  However,  being  able  to  show  convergence  even  up  to  a 
certain  propagation  distance  is  a  step  in  the  right  direction  toward  understanding 
the  behavior  of  these  collapses. 

Details  of  Collapse  Events 

With  time  convergence  of  the  collapse  events,  it  was  thought  that  the  size, 
shape,  and  behavior  of  these  occurrences  could  be  studied.  However,  convergence 
was  only  achieved  up  to  a  certain  propagation  distance.  This  distance  did  not 
allow  for  detailed  analysis  of  a  collapse’s  complete  evolution.  Although  it  is  clear 
that  the  collapses  continue  to  collapse  down  in  duration,  it  is  not  clear  if  this 
collapse  continues  indefinitely  or  stabilizes  at  a  certain  duration,  since  sampling 
frequency  could  not  accurately  describe  the  shape  and  size  after  a  certain  distance. 
Further  research  is  needed  to  understand  these  collapses  in  more  detail.  An  adap¬ 
tive  grid  that  can  have  more  resolution  where  the  collapse  events  are  occurring 
and  less  where  the  pulse  is  somewhat  stable  could  allow  for  this  type  of  analysis. 
The  FFT  approach  would  have  to  be  abandoned  in  the  time  dimension,  however, 
requiring  an  entirely  different  method  of  solving  the  propagation  equation  if  GVD 
is  included.  This  was  beyond  the  scope  of  this  research. 


53 


Stability  Analysis 


In  [21,23],  Niday  performed  an  analysis  on  the  stability  of  long  UV  pulses 
where  an  instantaneous  index  of  refraction  is  assumed.  He  studied  the  growth  of 
unstable  modes  with  respect  to  both  spatial  and  temporal  frequencies,  which  will 
be  referred  to  as  k±  and  fi,  respectively.  He  found  that  the  frequency  response  to 
the  instability  is  unbounded,  meaning  that  all  frequencies  present  in  the  numerical 
model  will  experience  growth,  and  this  growth  will  not  be  bandwidth  limited  as 
desired  for  a  spectral  code.  Thus,  he  concluded  that  long  pulses  suffer  from 
modulational  instability. 

Following  his  analysis,  we  will  discuss  the  stability  of  picosecond  pulses 
including  GVD  with  losses  considered  as  well.  We  start  with  the  propagation  and 
plasma  equations  (2.4)  and  (2.6),  including  only  the  imaginary  terms  (i.e.  no  loss 
terms)  and  a  simplified  plasma  equation,  giving 

d£  .1  2  .k"  d2£  .aujr  2  p41\ 

Tz  =,2k f  ■  !t w  -  l^rp£ + lkon^s  (  ’ 

%  =m2K  -  <*P2,  (4.2) 

where  b  =  We  then  assume  a  steady  state  solution  for  the  field  and  plasma. 
If  losses  were  included  at  this  point,  there  would  not  be  a  steady  state  solution 
because  the  pulse  would  continuously  decrease  due  to  the  losses.  The  steady  state 
solutions  are  solved  for  by  Niday,  who  shows  that  the  exact  solution  for  the  held 
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Figure  4.8:  Niday’s  steady  state  solution.  This  solution  contains  a  power  of 

roughly  470  MW,  and  beam  width  of  roughly  97  //rn.  The  peak  intensity  which 
defines  this  solution  is  3.2  x  1016  W/m2.  Also  shown  is  a  Gaussian  profile,  with 
the  same  peak  intensity  and  width.  From  reference  [23]. 

is  very  similar  to  a  Gaussian  (Figure  4.8).  Thus,  the  steady  state  solutions  for 
the  field  and  plasma  become 


£o  (r,z)  = 

=  u(r)eikz 

(4.3) 

Po(r)  = 

\T^\u(r)\K, 

V  a 

(4.4) 

respectively,  where  k  is  an  eigenvalue  that  can  be  solved  for  by  inserting  the  above 
equations  back  into  our  propagation  equation  (4.1).  This  results  in 


k 


aur 

~Y~ 


u(r)  |3  +  k0n2\u(r)\ 


2 


k0n3\u(r)\3  +  k0n2\u(r)  |2. 


(4.5) 
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Now,  if  we  let  £+,  £_,  and  p+  be  small  perturbations  compared  to  £  and  p, 


then  we  can  investigate  a  time-perturbed  solution  of 

S(r,z,t)  =  £0(r,z)  +  £+(r,z)e-m  +  £_(r,z)eint,  (4.6) 

p(r,  z,  t )  =  p0[r)  +  p+(r,  z)e~int  +  p+(r,  z)e*m,  (4.7) 


where  we  have  to  use  p+  and  p*+  since  the  plasma  density  is  real.  Following  Niday’s 
derivation  we  obtain 


Kb \£0\2K~2  (£*£+  +  £0£*_) 

Zctpo  —  iVt 

*  Kb\£0\2K~2  (£*0£^  +  £0£*+) 

2ap0  +  ^ 

Next,  we  insert  equations  (4. 6-4. 9)  into  equation  (4.1),  expand  the  terms,  and 
keep  only  the  linear  perturbation  terms.  In  this  step  we  must  utilize  the  fact  that 
K  =  3  for  UV  ionization  in  air.  We  then  get 

~7y  -  =  Vj_£o  +  ikoTi2 \£q \^£o  —  *  2  Po^o,  (4-10) 


(4.8) 

(4.9) 


for  the  zero  order  terms, 
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1  k" 
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■aUT  (  c  , 
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(4.11) 
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for  the  terms  oscillating  as  e  ,  and 


dS_ 

dz 


1  k" 

—  j_£_  +  i—£l~£_  +  ffco^2|^oh  [2£-  +  £+] 

z/c  z 


<tu;t  , 

—  (  Po£-  + 


Kb\£0 


2K 


2ap0  —  i£l 


[£-+£;]  (4.12) 


for  the  terms  oscillating  as  e  .  These  last  two  equations  give  the  evolution  of 
the  coupled  perturbation  fields  £+  and  £_.  We  now  assume  a  solution  for  these 
perturbation  fields  of 


£+(r,  z )  =  u. |_  exp(A  z  +  ikz  +  zkj_  •  r) 
£-(r,  z )  =  w_  exp(A*z  +  ikz  —  zkj_  •  r), 


(4.13) 


where  kj_  is  the  transverse  wave  vector  kxx  +  kyy ,  r  is  the  transverse  coordinate 
xx  +  yy ,  k  is  dehned  as  in  equation  (4.5),  and  A  now  represents  the  eigenvalue 
or  growth  rate  of  the  held.  Substituting  these  solutions  into  equations  (4.11)  and 
(4.12),  defining  the  contribution  from  the  as  k±  =  —  and  cancelling  the 
exponential  dependence  gives 


(ik  +  ik l  +  A)m+  =iGu+  —  iC  ( PqU+  +  AL(fl)  \u+  +  u*_]) 
+  iB  [2 u+  +  u*_]  , 

(ik  +  ik±  +  \*)u-  =iGu-  —  iC  (po«-  +  AL*(fl)  [u_  +  «+]) 
+  iB  [2 U-  +  , 


(4.14) 


(4.15) 
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where 


A(fi) 


Kb\£0\2K 
2apo  —  i£l 


C  = 


(JUT 

~Y~ 


B  =  k0n2\£o\2 
G(Q)  =  yfi2. 


(4.16) 


Putting  equation  (4.14)  and  the  complex  conjugate  of  (4.15)  into  matrix  form,  we 
have 
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u+ 

-Y 

-X 

u*_ 

ik  +  ikx_  +  A 


-ik  —  ik_ l  +  A 


u+ 


u 


(4.17) 


where 


X((2)  =  2 B  +  G{ ft)  -  C  (p0  +  A(fi)) 

(4.18) 

Y(Q)  =  B  —  CA(Q). 

(4.19) 

It  is  only  possible  to  have  a  nontrivial  solution  to  this  matrix  equation  if 


A2  =  Y2  -  (X  -k-  k±)2. 


(4.20) 


This  factors  as 


A 2  =  (B-  A(VL)Cf  -  (2 B  -  A{£l)C  +  G  -  k  -  k±  -  Cp0 )2  (4.21) 
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From  equation  (4.5)  and  equation  (4.16)  it  is  easy  to  infer 


k  =  B  —  Cpo . 


(4.22) 


Substituting  this  into  the  equation  for  A2  gives  a  reduced  form  of 


A2  =  (2 Y  -  k±  +  G )  (k±  -  G ) . 


(4.23) 


To  achieve  stability,  the  growth  rate,  Re( A),  needs  to  be  non-positive.  When 
GVD  is  not  included,  as  in  Niday’s  analysis,  this  only  occurs  at  zero  spatial 
frequency  or  zero  temporal  frequency.  This  is  shown  in  Figure  4.9  which  is  a 
reproduction  of  a  figure  from  Niday’s  dissertation  and  Figure  4.10  which  uses  the 
new  recombination  rate  and  the  higher  MPI  coefficient  of  3.9  x  10  33  m3/W2. 
Both  use  an  intensity  of  3.2  x  1016  W/m2.  With  these  two  figures  the  difference 
that  the  new  recombination  rate  makes  can  be  seen  in  the  maximum  value  of  the 
growth  rate  and  its  shape  in  the  k±  direction. 

When  GVD  is  included  very  little  changes  until  very  high  temporal  frequen¬ 
cies  are  reached.  Figures  4.9  and  4.10  were  reproduced  with  GVD  and  shown  at 
a  higher  temporal  frequency  in  Figures  4.11  and  4.12  so  the  difference  is  made 
clear.  Comparing  these  two  figures  shows  how  the  higher  MPI  coefficient  affects 
the  growth  rate.  The  higher  parameter  raises  the  value  of  the  growth  rate  at 
higher  spatial  frequencies.  It  should  be  noted  that  the  value  of  the  growth  rate 
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Figure  4.9:  Growth  rate  vs.  k±  and  fl  without  GVD  with  the  original  value  for 
a  and  (3^K\  for  an  intensity  of  3.2  x  1016  W/m2.  A  reproduction  of  Figure  4.2  in 
reference  [21]. 


Figure  4.10:  Growth  rate  vs.  k±  and  hi,  without  GVD,  with  an,  and  with  a 
higher  MPI  coefficient,  for  an  intensity  of  3.2  x  1016  W/m2. 
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(a)  Without  GVD 


(b)  With  GVD 

Figure  4.11:  (a)  Figure  4.9  shown  at  high  temporal  frequencies,  (b)  The  same 

figure  with  the  inclusion  of  GVD. 
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(a)  Without  GVD 


(b)  With  GVD 

Figure  4.12:  (a)  Figure  4.10  shown  at  high  temporal  frequencies,  (b)  The  same 

figure  with  the  inclusion  of  GVD. 


62 


is  not  zero  where  it  appears  to  be  negligible.  Figure  4.13  shows  that,  with  GVD 


included,  the  growth  rate  is  non-zero  at  all  non-zero  temporal  frequencies,  even 
when  k±  =  0. 

For  a  closer  look  at  the  stability  of  the  propagation  of  the  pulses,  let  us 
examine  the  high  frequency  limits  of  the  growth  rate.  As  k±  — >  oo , 


Re( A)  -»•  ±Im(Y). 


(4.24) 


Since  Y  is  not  dependent  on  G,  the  GVD  term,  this  is  the  same  result  that  Niday 
achieved  [21].  Thus,  there  will  always  be  a  positive  growth  rate,  although  it  will 
be  capped  for  large  spatial  frequencies.  This  limit  is  shown  in  Figure  4.14  for  the 
new  recombination  rate,  and  high  MPI  coefficient  showing  the  negative  root  as 
a  dashed  line.  However,  the  limit  of  infinite  spatial  frequency  is  never  reached 
because  of  the  paraxial  approximation  which  assumes  the  transverse  wave  vector, 
ki,  is  small  compared  to  kz.  Since  is  related  to  kj_,  its  maximum  value  is 
limited  to  approximately  3000  m-1  for  a  grid  with  Ax  =  0.8  mm. 

Taking  the  limit  of  high  temporal  frequencies,  as  D  — >  oo,  the  growth  rate 
remains  positive  without  an  upper  limit.  In  this  case,  Y  — »  B ,  so  with  GVD 
included, 


A2  ->  (2 B  —  k±  +  G)  (k l  -  G)  , 


(4.25) 
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Re(A)  (m_1) 


n(1012s_1) 


Re(A)  (m_1) 


(a)  Lower  Temporal  Frequencies 


n(1015s_1) 


Figure  4.13:  Growth  rate  vs.  when  k±  =  0  using  the  parameters  in  Table  3.1 
with  an  and  a  higher  MPI  coefficient,  (a)  Lower  temporal  frequencies,  (b)  Higher 
temporal  frequencies  where  A  =  248  nm  corresponds  to  7.6  x  1015  s-1 
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Figure  4.14:  Growth  rate  vs.  when  kj_  — >  oo  using  the  parameters  in 

Table  3.1  with  an.  The  dashed  line  is  the  negative  root  of  the  growth  rate  at 

this  limit. 

which  leads  to  an  instability  region  of  G(fl)  <  k±  <  2 B  +  G(f2)  where  the  growth 
rate  is  positive  even  as  11  — »  oo.  But  since  G  is  proportional  to  f22,  this  unstable 
region  continually  moves  to  higher  values  of  k±  and  decreases  in  size.  This  can 
be  seen  in  part  (b)  of  Figures  4.11  and  4.12.  When  not  including  GVD,  there  is 
a  constant  unstable  region  of  0  <  k±  <  2 B  as  determined  in  [21],  which  suggests 
that  difficulties  will  arise  when  using  numerical  techniques  for  analysis.  All  of 
the  frequencies  present  in  the  numerical  model  will  experience  growth,  and  the 
growing  held  will  not  be  bandwidth  limited,  which  is  desired  for  a  spectral  code. 
With  this  in  mind,  including  GVD  might  suggest  that  the  propagation  would  be 
more  stable  by  eliminating  the  growth  of  very  large  temporal  frequencies.  The 
results  of  our  simulations,  however,  reveal  breakdown  of  the  collapses  after  only 
short  propagation  distances  when  GVD  is  included  (shown  in  a  later  section). 
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This  may  be  due  to  the  fact  that  k±  =  0  does  not  lead  to  a  growth  rate  of  zero, 
but  when  GVD  is  neglected,  k±  —  0  does  lead  to  zero  growth. 

It  is  clear  that  in  all  scenarios,  pulse  propagation  is  inherently  unstable,  with 
a  positive  value  for  the  growth  rate  at  all  non-zero  frequencies.  Even  if  nonlinear 
losses  are  considered,  such  as  MPI  and  plasma  loss,  compared  to  the  linear  loss 
rate  of  2.5  x  10-4  m_1,  they  are  insignificant  at  low  intensities.  At  intensities 
corresponding  to  steady  state  solutions,  the  nonlinear  loss  rate  is  larger  than  the 
linear  loss  rate  at  0.045  m_1.  Only  when  intensities  approach  1018  W/m2  does 
the  nonlinear  loss  rate  compare  to  the  unstable  growth  rate. 


Propagation  Term  Values 

In  this  section,  a  closer  look  is  taken  at  each  of  the  individual  terms  in  the 
propagation  equation  (2.4)  where 


d£  _  ik"  d2£ 
~dz~~  ~YW 

^L\£]2K-2£ 

~  |  (!  +  ivr)  p£ 
ik0n2\£\2£ 


(4.26) 


will  be  referred  to  as  the  GVD,  MPI,  plasma,  and  Kerr  terms,  respectively.  Note 
that  the  diffraction  term  is  not  listed  in  this  discussion  as  it  applies  more  to 
the  spatial  dimensions  versus  the  temporal  dimension.  Using  the  propagation 
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code  and  Mathematica®,  the  values  of  these  terms  were  determined  at  specific 
temporal  points  along  the  propagation  axis  (x=0,  y= 0,  t).  To  understand  which 
terms  played  a  major  role  in  forming  the  collapse  events,  the  coordinate  of  the 
maximum  value  of  the  intensity  was  found.  Then,  the  value  of  the  terms  was 
determined  slightly  before,  at,  and  slightly  after  this  coordinate.  A  run  with 
the  new  recombination  coefficient,  Nt  =  16384,  and  At  =  0.9  fs  was  used  in 
determining  these  values.  A  time  slice  where  x  =  0  of  all  of  these  terms  is  shown 
in  Figures  4.15  through  4.17  at  various  propagation  steps.  The  units  for  these 
terms  are  the  normalized  field  units  per  meter.  At  the  beginning  of  propagation, 
z  =  2.5  cm,  the  first  peak  has  barely  started  to  form  at  t  —  —3.265  ps.  On  either 
side  of  this  peak,  t  =  —3.355  ps  and  t  =  —3.175  ps,  there  is  little  difference  in  the 
terms.  The  largest,  the  plasma  term,  at  0.101  for  the  real  part  increased  to  0.117, 
and  the  smallest  term,  the  Kerr  term,  at  -10.5123  for  the  real  part  increased  to 
-10.5015.  The  real  part  of  the  MPI  term  was  negative  and  an  order  of  magnitude 
smaller  than  the  plasma  term,  and  the  real  part  of  the  GVD  term  was  three  orders 
of  magnitude  smaller  than  the  plasma.  The  summation  of  the  Kerr  and  plasma 
terms  did  increase  by  0.02  after  the  peak. 

However,  after  the  collapse  has  had  a  chance  to  form  at  a  propagation 
distance  of  15  cm,  the  behavior  of  the  terms  is  noticeably  different  in  front  of  and 
behind  the  peak.  The  peak  is  now  at  t  —  —4.5  ps.  Just  slightly  before  the  peak 
at  t  —  —4.52  ps  the  real  part  of  the  plasma  term  is  —74.43,  and  slightly  after  the 
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Figure  4.15:  Time  slice  of  the  terms  in  the  propagation  equation  at  2.5  cm. 

Both  a  full  view  (a)  and  a  magnification  (b)  of  the  absolute  values  of  each  term 
are  shown  for  a  run  with  the  new  recombination  coefficient  and  high  resolution. 
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Figure  4.16:  Time  slice  of  the  terms  in  the  propagation  equation  12.5  cm.  Both 
a  full  view  (a)  and  a  magnification  (b)  of  the  absolute  values  of  each  term  are 
shown  for  a  run  with  the  new  recombination  coefficient  and  high  resolution. 
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Figure  4.17:  Time  slice  of  the  terms  in  the  propagation  equation  at  50  cm. 

Both  a  full  view  (a)  and  a  magnification  (b)  of  the  absolute  values  of  each  term 
are  shown  for  a  run  with  the  new  recombination  coefficient  and  high  resolution. 
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peak  at  t  —  —4.48  ps  it  is  —1880.  However,  the  real  part  of  the  Kerr  term  only 
increases  from  1170  to  1730  making  the  summation  of  the  two  highly  positive 
before  the  peak  and  negative  after  the  peak.  This  just  perpetuates  the  increasing 
slope  on  the  front  of  the  collapse  and  decreasing  slope  on  the  back,  causing  the 
collapse  event  to  continue  to  compress  temporally.  The  GVD  and  MP1  terms 
still  remain  orders  of  magnitude  smaller  than  the  other  terms,  thus  they  do  not 
contribute  much  to  this  analysis. 

In  order  to  abate  this  high  slope  on  the  front  of  the  pulse,  the  plasma  would 
need  to  be  higher  at  the  beginning  of  the  pulse.  This,  however  is  a  difficult  task 
to  achieve  when  working  with  just  one  pulse,  since  the  plasma  is  created  after  the 
leading  edge  of  the  pulse  has  passed.  This  is  why  studies  dealing  with  multiple 
pulses  is  desired.  It  is  suggested  that  a  leading  pulse  can  create  a  plasma  for  the 
following  pulse  to  enter,  thus,  there  is  already  an  established  plasma  term  at  the 
beginning  of  the  pulse  that  would  be  able  to  balance  out  the  Kerr  term  creating 
a  stable  pulse  as  opposed  to  a  collapse. 

As  discussed  in  the  section  on  GVD  in  Chapter  II,  group  velocity  dispersion 
was  included  in  this  study  of  the  collapse  events  of  UV  pulses  due  to  the  extremely 
short  time  scale  of  the  collapses.  It  was  thought  that  GVD  might  be  able  to 
spread  these  collapses  out,  thus  curbing  their  growth  and  compression.  Although 
adding  the  GVD  term  did  affect  the  collapse  peaks  slightly,  it  also  caused  extreme 
perturbations  of  the  collapses.  The  nature  of  the  second  derivative  of  the  field 
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with  respect  to  time  found  in  the  GVD  term,  is  to  magnify  any  noise  in  the  held. 
A  plot  of  the  absolute  value  of  the  GVD  term  for  z  =  2.5  cm  is  shown  in  Figure 
4.18  revealing  the  magnification  of  noise.  The  GVD  parameter,  k"  was  varied  to 
help  determine  how  much  the  GVD  term  affected  the  pulse  propagation  code.  It 
was  found  that  the  higher  the  value  of  k" ,  the  sooner  the  collapsed  peaks  showed 
signs  of  perturbation.  This  is  shown  in  Figure  4.19.  It  displays  a  magnification 
of  a  collapse  event  of  a  pulse  having  the  same  parameters  from  Table  3.1,  except 
for  k"  being  0,  1.45  x  1CT28,  and  1.45  x  10“26  s2/m  with  Nt  =  4096,  At  =  3  fs. 
The  old  value  for  a  is  used.  Without  GVD  included,  the  pulse  can  propagate  to  a 
distance  of  75  cm  and  farther,  but  when  normal  GVD  is  included,  perturbations 
appear  at  47.5  cm.  For  the  bigger  GVD  value,  the  pulse  only  gets  to  25  cm  before 
it  shows  signs  of  perturbations.  To  avoid  these  perturbations  in  the  results  of  the 
code  due  to  numerical  error  being  magnified,  we  neglected  GVD  for  many  of  the 
simulations. 

The  terms  in  the  plasma  propagation  equation  (2.5)  were  also  investigated. 
The  diffusion  strength  D  is  still  assumed  to  be  zero,  but  the  avalanche  ionization 
term  is  included  to  check  the  validity  of  neglecting  it  in  the  main  propagation 
code.  The  following  terms  will  be  referred  to  as  the  MPI,  recombination  and 
avalanche  terms,  respectively: 


Cp\£\\ 


f3W\£\2K 

Khw 


ap2, 


(4.27) 
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Figure  4.18:  A  middle  slice  of  the  absolute  value  of  the  GVD  term  with  respect 
to  time.  Shown  for  the  pulse  at  z  =  2.5  cm  with  the  new  alpha  value.  Note  the 
perturbations  due  to  taking  the  second  derivative  numerically. 

where  C  —  a/ Eg  (assuming  a  background  index  of  one)  and  the  ionization  energy 
Eg  ps  11  eV  [22],  As  shown  in  Figure  4.20  the  MPI  term  is  much  greater  than  the 
avalanche  term,  so  this  assumption  stands. 


Experimental  Verification 

In  reference  [24],  Chalus  et  al.  report  experimentally  confirming  the  forma¬ 
tion  of  a  filament  with  the  following  conditions:  power  of  500  MW,  beam  waist 
of  100  /jin,  pulse  duration  of  200  ps.  These  parameters  were  used  as  initial  con¬ 
ditions  for  the  propagation  code  with  both  a  supergaussian  time  profile  and  a 
regular  Gaussian  profile.  The  supergaussian  temporal  shape  of  the  pulse  is  cho¬ 
sen  to  simulate  the  time  independent  solution  assumed  by  the  theory  proposed 
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Figure  4.19:  Part  (a)  shows  the  first  collapse  of  the  pulse  with  no  GVD  term  at 
z  =  0.75  m.  Part(b)  shows  the  first  collapse  of  the  pulse  with  k"  =  1.45  x  10~28 
s2/m  for  the  GVD  term  at  z  =  0.475  m.  Part(c)  shows  the  second  collapse  of  the 
pulse  with  k"  =  1.45  x  10~26  s2/m  for  the  GVD  term  at  z  —  0.275  m.  The  first 
collapse  was  already  perturbed  to  such  high  frequency,  aliasing  was  occurring. 
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Figure  4.20:  A  middle  slice  of  the  plasma  propagation  terms  with  respect  to 

time.  Shown  for  the  pulse  at  z  —  2.5  cm  with  both  the  new  and  old  alpha  values. 

to  explain  long  distance  propagation  [11,21,23].  The  new  recombination  rate  was 
used  as  well  as  the  higher  MPI  coefficient  to  ensure  convergence.  The  superguas- 
sian  run  resulted  in  similar  collapse  events  as  previously  discussed.  However,  the 
Gaussian  run  resulted  in  the  pulse  somewhat  collapsing,  but  into  only  one  event 
that  stabilized  over  the  first  70  cm  of  propagation.  It  could  be  able  to  propa¬ 
gate  stably  for  even  longer  distances.  This  could  not  be  determined  due  to  time 
constraints,  however.  Figures  4.21  through  4.25  show  a  3D  representation  of  this 
Gaussian  pulse  as  it  propagates,  along  with  its  /e-space  form  in  (a)  as  well  as 
the  fluence  profile  of  the  pulse  in  (b).  Figure  4.26  shows  a  time  slice  along  the 
x-axis  as  well  as  the  fluence  profile  of  the  same  pulse.  A  larger  spatial  grid  size 
of  Nx  =  Ny  =  512  has  prevented  the  fluence  profile  from  showing  signs  of  grid 
overflow.  Note,  however,  that  by  50  cm  the  pulse  is  starting  to  diverge  beyond 
the  limits  of  the  spatial  grid.  The  ABC’s  will  not  be  able  to  absorb  the  held  as  it 
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diverges  further  off  the  grid  at  farther  propagation  distances.  When  the  fluence 
profile  begins  to  show  signs  of  this,  the  results  from  the  propagation  code  become 

invalid. 

Summary 

In  summary,  time  convergence  is  achieved  up  to  at  least  1  m  of  propagation 
when  a  larger  value  of  /3 is  used.  However,  with  the  original  value,  convergence 
can  only  be  achieved  up  to  25  cm  even  with  Nt  =  16384  and  Af=0.9  fs.  This 
prevents  a  detailed  analysis  of  the  collapse  events.  Only  the  fact  that  the  collapses 
increase  in  intensity  and  decrease  in  duration  during  that  first  25  cm  can  be  con¬ 
firmed.  In  order  to  study  this  further,  one  would  need  a  code  that  has  adaptable 
grid  parameters  to  achieve  enough  resolution  of  the  collapses,  yet  not  overload  the 
computer’s  ability,  or  one  would  need  access  to  a  supercomputer.  It  was  shown 
using  a  stability  analysis  that  in  all  cases  the  pulse  remains  modulationally  unsta¬ 
ble.  The  GVD  term  specifically  introduces  instability  for  all  temporal  frequencies 
at  k±  =  0.  The  high  temporal  frequency  limit  was  taken  to  reveal  the  effect 
that  GVD  has  on  the  numerical  stability,  and  the  high  spatial  frequency  limit 
revealed  the  same  information  as  in  reference  [21].  The  second  derivative  in  the 
GVD  term  is  believed  to  cause  nonphysical  perturbations  to  show  in  the  pulse 
and  collapses  after  just  a  short  propagation  distance  during  runs  including  GVD. 
A  closer  look  at  the  actual  values  of  the  propagation  terms  revealed  that  the  Kerr 
and  plasma  terms  interacted  to  determine  the  slope  of  the  field.  The  other  terms 
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z=1.25  cm 


Figure  4.21:  The  propagation  of  a  Gaussian  200  ps  duration  pulse  along  with 
its  fc-space  representation.  (b)is  the  fluence  profile.  The  peak  amplitude  of  (a) 
was  normalized  to  one  at  z  =  0,  and  the  following  plots  are  on  the  same  scale. 
The  fluence  profiles  in  the  following  plots  are  not  to  scale. 
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z=15  cm 


Figure  4.22:  The  propagation  of  a  Gaussian  200  ps  duration  pulse  along  with 
its  fc-space  representation.  (b)is  the  fluence  profile. 
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z  =  3  0  cm 


(b) 


Figure  4.23:  The  propagation  of  a  Gaussian  200  ps  duration  pulse  along  with 
its  fc-space  representation.  (b)is  the  fluence  profile. 
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z=50  cm 


z=50  cm 


(b) 


Figure  4.24:  The  propagation  of  a  Gaussian  200  ps  duration  pulse  along  with 
its  /c-space  representation.  (b)is  the  fluence  profile.  Note  that  the  pulse  in  (a)  is 
starting  to  spread  beyond  the  limits  of  the  spatial  grid.  It  has  not  effected  the 
fluence  profile,  yet. 
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z=70  cm 
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z  =  7  0  cm 


(b) 

Figure  4.25:  The  propagation  of  a  Gaussian  200  ps  duration  pulse  along  with 
its  fc-space  representation.  (b)is  the  fluence  profile.  Note  that  the  pulse  in  (a) 
is  continuing  to  spread  beyond  the  limits  of  the  spatial  grid.  It  still  has  not 
effected  the  fluence  profile.  Further  propagation  steps  were  unavailable  due  to 
time  constraints. 
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Figure  4.26:  A  time  slice  of  a  Gaussian  200  ps  duration  pulse  at  (a)  1.25  cm, 
(b)  15  cm,  (c)  30  cm,  (d)  50  cm,  and  (e)  70  m. 
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were  orders  of  magnitude  smaller  and  therefore  contributed  little  to  the  evolu¬ 
tion  of  the  collapses.  It  is  suggested  that  if  a  plasma  could  be  introduced  to  the 
front  of  the  pulse,  say  by  a  previous  pulse,  the  plasma  term  could  balance  out  the 
Kerr  term  and  establish  a  stable  pulse  as  opposed  to  collapse  events.  The  terms 
in  the  plasma  propagation  were  also  examined  to  confirm  the  approximation  of 
neglecting  avalanche  ionization. 


83 


V.  Conclusions  and  Suggestions 


This  chapter  will  discuss  the  conclusions  that  can  be  drawn  from  the  research 
presented  and  explains  the  impact  this  research  has  on  the  scientific  investigation 
of  hlamentation  and  future  applications  for  the  United  States  Air  Force.  It  also 
gives  suggestions  for  future  research  that  can  stem  off  what  has  been  done  here. 

Research  Conclusions 

In  addition  to  Niday’s  research  on  the  modeling  of  the  propagation  of  long 
UV  pulses,  Muller  came  to  several  relevant  conclusions.  First,  he  established 
numerical  grid  parameters  in  the  spatial  dimension  that  accurately  displayed  the 
pulse’s  propagation.  He  could  not,  however,  establish  grid  parameters  for  the 
time  dimension  that  prevented  numerical  sampling  errors,  because  the  collapses 
condensed  to  durations  on  the  order  of  the  smallest  grid  spacing  he  was  able  to 
attain.  Second,  he  established  that  the  collapses  are  strongly  related  to  the  plasma 
generated  by  the  pulse.  Third,  he  discovered  that  the  input  power  of  the  pulse 
was  directly  proportional  to  the  rate  at  which  the  collapses  formed.  Working  from 
this  information  and  the  foundation  laid  by  Niday,  this  work  presents  several  more 
conclusions  on  the  phenomenon  of  collapse  events  forming  on  the  leading  edge  of 
a  self-focusing  UV  picosecond  pulse. 

The  first  conclusion  of  this  research  is  that  the  collapse  events  become  too 
short  in  duration  for  the  computational  power  available  to  be  able  to  handle  the 
numerical  parameters  needed  to  achieve  convergence  of  these  events.  Convergence 
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was  only  demonstrated  for  the  first  25  cm  of  propagation  before  the  collapses 
condensed  to  a  duration  unable  to  be  represented  by  a  grid  spacing  of  0.9  fs. 
This  initial  convergence  confirms  that  at  least  the  formation  of  the  collapses  is 
being  appropriately  represented  by  the  numerical  code,  even  if  their  complete 
evolution  remains  a  mystery.  Complete  time  convergence  was  achieved  up  to  1  m, 
however,  using  a  MPI  coefficient  10  times  larger  than  originally  specified.  This 
larger  value  causes  the  losses  to  overcome  the  increasing  intensity  of  the  collapse 
so  the  numerical  parameters  are  able  to  represent  it  longer  before  it  eventually 
compresses  to  a  duration  on  the  order  of  our  grid  spacing.  This  confirms  that 
the  collapses  continue  to  decrease  in  duration,  but  do  not  continually  increase  in 
intensity,  at  least  with  this  higher  MPI  loss. 

The  second  conclusion  established  in  this  research  stems  from  the  investiga¬ 
tion  of  the  individual  propagation  term  values.  The  ever  increasing  front  slope  of 
the  collapse  is  due  to  the  Kerr  term  dominating  the  plasma  term.  The  decreas¬ 
ing  negative  slope  on  the  trailing  edge  of  the  collapse  is  due  to  the  plasma  term 
decreasing  greatly  while  the  Kerr  term  only  increases  slightly.  The  other  terms  in 
the  propagation  equation  are  orders  of  magnitude  smaller  than  these  two  terms. 
This  interaction  suggests  that  if  more  plasma  is  present  at  the  front  of  the  pulse, 
the  plasma  term  will  be  more  negative,  possibly  allowing  for  a  more  neutral  slope 
of  the  leading  edge  of  the  pulse,  establishing  stable  propagation  of  the  pulse  as 
opposed  to  the  collapse  events  being  generated. 
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Third,  it  was  concluded  that  GVD,  although  affecting  the  collapses  slightly 
due  to  their  short  duration,  can  be  neglected.  The  nature  of  the  GVD  term  induces 
an  unbounded  growth  of  all  frequencies  introduced  by  the  numerical  code.  By 
including  GVD  in  the  analysis,  any  and  all  noise  is  amplified  causing  nonphysical 
perturbations  in  the  results  for  propagation.  The  fact  that  the  value  of  the  GVD 
term  in  the  propagation  term  investigation  was  insignificant  compared  to  the  Kerr 
and  plasma  terms  and  the  results  of  the  stability  analysis  lead  to  this  conclusion. 

Another  conclusion  of  this  study  also  stems  from  the  stability  analysis.  It 
was  discovered  that  linear  instability  is  inherent  in  the  propagation  of  UV  pulses 
with  the  initial  conditions  specified  in  this  study.  Even  zero  spatial  frequencies 
have  a  positive  growth  rate  when  GVD  is  included.  Losses  were  found  to  compare 
to  this  growth  rate  only  at  high  intensities  on  the  order  of  1018  W/m2.  This  is 
found  to  be  true  with  all  combinations  of  the  original  and  new  values  for  the 
recombination  rate  and  the  MPI  coefficient. 

Finally,  the  initial  conditions  of  an  experiment  that  generated  hlamentation 
was  introduced  to  the  propagation  code  using  the  new  recombination  rate  and  a 
higher  MPI  coefficient.  This  run  generated  only  one  collapse  event  that  stabilized 
and  continued  to  propagate  farther  than  1  m.  However,  results  after  30  cm  start 
to  exhibit  signs  that  at  least  part  of  the  pulse  has  spread  beyond  the  limits  of  the 
spatial  grid,  introducing  numerically  induced  phenomena. 


Impact  of  Research 


The  study  of  the  filamentation  of  UV  light  is  extremely  relevant  to  the 
United  States  Air  Force.  If  we  can  understand  the  intricate  details  behind  how 
light  filaments  are  able  to  travel  long  distance  without  excessive  dispersion,  diffrac¬ 
tion  or  even  annihilation  due  to  obscurants,  our  ability  to  harness  this  phenomenon 
for  numerous  applications  is  greatly  increased.  By  modeling  filamentation  accu¬ 
rately,  the  ideal  initial  conditions  can  be  determined  that  promote  desired  effects 
for  applications  such  as  remote  sensing  of  aerosols  in  the  atmosphere,  protection 
of  equipment  and  personnel  from  lighting  strikes,  or  the  ability  to  direct  energy  at 
certain  targets  for  their  destruction.  Niday’s  UV  propagation  model  is  a  piece  of 
this  understanding.  Using  it  to  investigate  the  behavior  of  how  light  travels  and 
the  nonlinear  effects  that  counteract  normal  diffraction,  we  are  one  step  closer  to 
the  realization  of  these  beneficial  applications  to  the  Air  Force. 

By  achieving  time  convergence  of  the  pulse  propagation  up  to  a  certain 
distance,  we  are  able  to  confirm  that  the  results  presented  here  are  not  due  to 
numerical  error,  but  accurately  represent  the  behavior  of  the  model.  However, 
until  the  model  accurately  represents  experimental  data,  we  must  search  for  the 
right  initial  conditions  or  a  better  model.  The  numerical  error,  found  by  previ¬ 
ous  research,  of  collapse  events  occurring  rapidly  on  the  trailing  edge  after  some 
distance  of  propagation  has  been  eliminated.  The  benefit  of  this  convergence  and 
elimination  of  numerical  error  is  validation  of  further  results  involving  the  intrica- 
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cies  of  propagation.  The  failure  of  convergence  at  further  distances  also  points  out 
the  need  for  an  adaptable  grid  size  or  supercomputer  that  can  achieve  even  more 
resolution  of  the  collapses  as  they  continue  to  intensify  and  shorten  in  duration. 

Having  a  more  in  depth  analysis  of  the  GVD  parameter  provides  several 
benefits  to  the  continued  analysis  of  UV  propagation.  The  stability  analysis  re¬ 
vealed  that  including  GVD  brought  stability  in  the  high  temporal  range,  but  also 
brought  instability  in  the  low  spatial  range,  which  could  be  the  cause  of  additional 
perturbations  in  the  pulse  during  propagation.  By  combining  the  understanding 
of  why  fragmentation  of  the  collapses  occurs  in  relation  to  GVD  and  the  real¬ 
ization  that  GVD  does  not  play  a  big  enough  role  to  curb  the  time  shortening 
of  the  collapse  events,  we  can  neglect  GVD  in  future  research  dealing  with  rel¬ 
atively  long  pulses.  This  simplifies  the  propagation  equation  for  easier  analysis 
and  manipulation. 

The  results  of  the  investigation  on  the  values  of  each  of  the  propagation 
terms  allow  future  research  to  focus  on  the  terms  which  play  the  biggest  roles  in 
the  formation  of  these  collapse  events.  The  interaction  of  these  terms,  Kerr  and 
plasma,  determined  the  sign  of  the  slope  of  the  field.  Knowing  this  suggests  that 
having  a  preceding  pulse  that  establishes  a  higher  plasma  term  at  the  front  of 
the  pulse  could  lead  to  the  elimination  of  the  collapse  events  that  occur  on  the 


leading  edge  of  the  pulse. 


The  stability  analysis  done  in  this  research  shows  that  with  GVD  considered 
and  all  combinations  of  recombination  rates  and  MPI  coefficients,  there  is  still 
unstable  growth  for  all  spatial  and  temporal  frequencies.  This  is  valuable  in  that 
it  raises  the  question  of  whether  we  have  correct  values  for  all  of  our  parameters 
or  if  hlamentation  is  even  theoretically  possible  with  the  model  assumed  for  these 
simulations.  Possibly,  hlamentation  could  arise  from  the  instability  itself.  With 
the  knowledge  of  this  instability  we  are  one  step  closer  to  finding  the  answers  to 
the  filament  phenomenon. 

With  the  simulation  of  the  propagation  of  a  pulse  imitating  one  that  ex¬ 
perimentally  created  a  filament,  the  code  produced  only  one  collapse  event  that 
was  able  to  converge.  The  results  of  this  run  could  be  consistent  with  the  ex¬ 
perimental  results,  thereby  validating  the  propagation  code.  If  this  is  indeed  the 
case,  attention  can  turn  toward  the  conditions  needed  to  produce  these  results 
as  opposed  to  determining  the  details  of  the  non-successful  runs.  However,  the 
numerical  grid  size  still  limited  these  results. 

Suggestions  for  Further  Research 

Because  of  the  time  limitation  on  this  research,  there  are  many  areas  in 
the  long  range  propagation  of  UV  laser  pulses  that  have  been  left  unaddressed. 
This  research  has  also  raised  many  new  questions  that  should  be  answered  about 
the  interaction  of  self-focusing,  diffraction,  and  other  losses.  Further  expansion 
of  Niday’s  propagation  model  is  needed  to  be  able  to  fully  grasp  the  concepts  of 


UV  filamentation  and  to  better  represent  the  unknown  intricacies  of  leading  edge 
collapse  events.  Experimental  verification  is  also  needed  to  be  done  alongside  this 
theoretical  modeling  work  in  order  to  confirm  its  results. 

Expansions  to  the  current  propagation  model  could  involve  inclusion  of  cer¬ 
tain  terms  that  were  neglected  for  simplification.  One  of  these  terms  is  the  last 
term  in  equation  (2.2),  the  Raman  effect.  Including  this  scattering  turns  n2  into 
a  non- instantaneous  self-focusing  term  by  introducing  a  time  dependence  part 
according  to 


ik0n2\S\2S  — >  ik0(  1  -  f)n2\S\ 2£  +  ik0fn2 
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where  /  is  the  fraction  of  the  Kerr  effect  which  is  not  instantaneous,  but  delayed 
due  to  stimulated  molecular  Raman  scattering.  Diffusion,  D  in  the  plasma  equa¬ 
tion  (2.5),  could  also  be  included  by  making  it  nonzero.  Another  way  to  expand 
the  code  for  better  propagation  is  to  add  Rayleigh  scattering  once  greater  dis¬ 
tances  are  achieved.  This  would  attenuate  the  power  of  the  beam  with  respect  to 
propagation  distance  by 


dP 

dz 


_p(K)l_p3  _  n^4  J_P5/ 2  _  aRp 
w 4  w6 


(5.2) 


where  the  first  term  on  the  right  hand  side  is  the  power  attenuation  through  the 
photon  ionization,  the  second  represents  plasma  absorption,  and  the  third  is  the 
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power  loss  through  Rayleigh  scattering.  For  a  wavelength  of  248  nm,  olr  =  5  x  10-4 
m_1  [17].  If  propagating  distances  of  kilometers,  one  could  also  add  variation  to 
the  background  index  of  refraction  no,  introducing  dependence  on  turbulence, 
altitude,  and  so  forth. 

Another  improvement  to  the  model  code  that  would  aid  in  the  study  of  the 
leading  edge  collapse  events  is  to  introduce  a  variable  time  step.  By  being  able  to 
use  a  larger  At  where  the  pulse  is  fairly  stable,  but  at  the  same  time  use  a  very 
small  time  step  only  where  the  collapse  events  are  occurring,  one  could  keep  the 
grid  size  at  a  reasonable  level  while  still  achieving  high  resolution  in  the  critical 
areas.  It  could  be  implemented  by  basing  the  time  step  on  the  derivative  of  the 
held  itself  instead  of  being  a  constant  input  by  the  user.  If  the  derivative  of  the 
held  is  small,  At  can  be  larger,  but  where  the  derivative  is  large,  the  time  step 
would  be  decreased  allowing  for  more  resolution  at  the  collapses.  Further  thought 
is  needed  to  consider  how  the  grid  size  would  be  determined  for  this  variable  time 
step.  This  adaptation  to  the  code  would  allow  for  much  better  resolution  while 
requiring  less  time  and  memory  to  compute.  This  higher  resolution  could  lead  to 
more  insight  as  to  the  nature  of  these  collapse  events.  However,  the  FFT  could 
no  longer  be  taken  in  the  time  dimension  since  it  requires  a  hxed  time  step.  This 
is  possible  if  GVD  is  not  considered  in  the  equation.  Another  option  to  increase 
resolution,  but  also  avoid  the  numerical  error  introduced  by  the  limited  grid  sizes, 
is  to  run  the  code  on  supercomputers.  This  could  be  the  key  to  bridging  the  gap 
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between  experimental  results  and  the  results  produced  numerically  with  the  3D 
propagation  code. 

More  detailed  study  of  each  collapse  event  could  also  be  done  by  taking  a 
collapsed  section  of  a  pulse  that  has  been  propagated  by  the  original  code  and 
using  it  as  an  input  to  a  higher  resolution  propagation  code.  The  higher  resolution 
would  have  to  be  interpolated  from  the  less  frequent  data  points  of  the  original  run. 
Care  would  have  to  be  taken  to  extract  the  collapse  event  from  the  original  run 
before  it  has  collapsed  to  a  size  that  the  grid  can  no  longer  accurately  represent  so 
that  the  interpolated  higher  resolution  points  will  be  correctly  determined.  With 
this  idea,  one  could  see  if  these  individual  collapse  events  are  able  to  propagate 
by  themselves  after  splitting  from  the  pulse. 

The  introduction  of  obscurants  to  the  pulse  path  could  be  another  extension 
of  this  propagation  model.  By  zeroing  out  a  portion  of  the  field  at  a  particular 
propagation  step,  an  obscurant  can  be  simulated.  Studies  could  be  done  to  see 
how  certain  pulses  overcome  or  are  dissipated  by  this  obstruction.  Experiments 
could  be  done  on  different  sizes  and  transparencies  of  the  simulated  barriers,  as 
well  as  how  placement  affects  the  outcome.  Does  placing  the  obstruction  on  axis 
or  off,  or  at  the  beginning  of  propagation  or  after  collapse  events  have  already 
been  allowed  to  form  make  a  difference  in  the  pulse  being  able  to  continue  to 
propagate?  All  of  these  scenarios  would  help  in  understanding  the  self-healing 
properties  of  filaments. 


92 


Because  the  plasma  density  is  not  immediately  terminated  at  the  end  of  the 
pulse,  but  dies  down  slowly,  the  idea  of  introducing  multiple  pulses  to  the  code  is 
intriguing.  Would  there  be  enough  left  over  plasma  that  would  interact  with  the 
next  laser  pulse  causing  a  completely  different  effect  than  was  seen  with  the  first 
pulse?  Maybe  the  plasma  would  just  cause  an  increase  or  decrease  in  the  effects 
presented  in  this  study.  After  many  pulses  have  been  propagated,  the  interaction 
of  them  could  lead  to  steadier  reaction  of  the  following  pulses  eliminating  the 
leading  edge  collapses.  Enough  questions  are  raised  when  thinking  about  pulse 
interactions  to  warrant  an  investigation. 

The  idea  that  X-waves  could  be  a  representation  of  filaments  propagating 
long  distances  is  an  exciting  and  promising  idea.  Research  could  be  done  to 
investigate  how  X-waves  could  relate  to  these  collapse  events  found  on  the  leading 
edge  of  pulses  being  propagated  with  the  model  presented  here.  Studies  could  be 
done  on  what  conditions  would  be  necessary  to  allow  these  collapses  to  form 
into  X-waves.  An  X-wave  could  also  be  used  as  an  initial  input  to  study  how 
the  simulation  allows  it  to  propagate.  Once  a  true  X-wave  is  able  to  propagate 
without  undergoing  diffraction  or  dispersion,  slight  variations  could  be  made  to 
the  input  pulse  to  deviate  from  the  ideal  conditions  to  see  if  the  deviations  could 
be  overcome  as  the  pulse  propagates  to  form  this  X-wave. 

More  experimental  data  is  needed  with  regard  to  longer  UV  pulses  at  high 
powers.  The  successes  and  failures  of  these  experiments  could  help  guide  following 
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research  to  the  areas  that  are  most  feasible  and  applicable  to  UV  radiation.  An 
extension  of  the  stability  analysis  that  includes  nonlinear  terms  could  also  be 
conducted  to  determine  if  stability  is  established  when  higher  order  terms  are 
considered.  Also,  a  more  in  depth  loss  analysis  could  be  done. 

Further  studies  could  be  done  with  regard  to  the  needs  of  the  Air  Force  in  the 
area  of  UV  long  distance  propagation.  A  survey  could  be  done  to  determine  the 
specific  areas  that  deal  with  remote  sensing,  laser  induced  and  guided  lightning,  or 
directed  energy  applications  in  which  the  military  is  particularly  interested.  From 
this,  one  could  receive  guidelines  for  wavelengths,  powers,  and  other  specifications 
required  for  use  in  combat  to  use  as  initial  conditions  in  simulations. 

Finally,  more  runs  should  be  made  with  the  initial  conditions  described 
in  the  experiment  of  reference  [24]  to  include  larger  spatial  grid  sizes  or  higher 
ABC’s  to  alleviate  the  restriction  of  the  pulse  due  to  numerical  limitations.  Also, 
the  wavelength  could  be  changed  from  248  nm  to  266  nm  as  was  used  in  the 
experiment.  If  these  runs  are  indeed  successful  in  creating  a  filament  that  can 
propagate  long  distances,  by  starting  with  this  scenario,  small  changes  can  be 
made  to  pinpoint  exactly  what  is  causing  the  code  to  exhibit  the  multiple  non- 
convergent  collapses  whether  it  be  grid  size,  values  of  the  recombination  rate  or 
MP1  coefficients,  or  the  initial  shape  of  the  pulse. 
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